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NOMENCLATURE 


constants used to represent the r distribution in the downwash 
sect ion. 

vector from a vortex filament to a point in space. 

wing span in downwash section or wake width in boundary layer 
sect i on . 

wing drag coefficient, 
wing lift coefficient, 
wing chord. 

wing section drag coefficient, 
wing section lift coefficient. 

eddy viscosity parameter in boundary layer solution, 
total drag. 

cylinder diameter in boundary layer section, 
profile drag per unit wing span. 

function of £ and n which will appear as the unknown in the 
transformed boundary layer equations. 

turbulent boundary layer form factor. 

cylinder length in' boundary layer section or perpendicular distance 
from the point of interest to vortex filament in downwash section. 

vector having direction of infinitesimal vortex filament. 

streamwise index in boundary layer finite difference solution. 

wing section lift curve slope. 

normal index in boundary layer finite difference solution, 
dynamic pressure. 

unknown in finite difference boundary layer recursive relation, 
wing thickness. 



NOMENCLATURE (Continued) 


U - frees+ream velocity in boundary layer section. 

u - induced freestream velocity component in the downwash section or 

velocity component parallel to streamline direction in boundary 
layer section. 

V - freestream velocity vector in downwash section. 

v - induced spanwise velocity vector In downwash section or normal 

velocity. 

W - total velocity induced at a point by a vortex. 

w - induced total velocity vector in general downwash section or induced 

vertical velocity component. 

x - freestream coordinate in downwash section or streamline coordinate 

in boundary layer section.- 

y - spanwise coordinate in the downwash section or normal coordinate in 

boundary layer section. 

z - vertical coordinate in the downwash section. 

a - angle of attack. 

r - strength of the bound vortex at any wing section. 

e - downwash angle in downwash section or eddy viscosity in boundary 

layer section. 

z, - nondimensiona I' wake width. 

r) - nond imens Iona I normal coordinate in boundary layer section or 

normalized loss of wake dynamic pressure. 

A - sweep angle. 

£ - nondimensiona I distance behind wing trailing edge. 

p - dens i ty of air. 

r - turbulent shearing stress. 

ip - two-dimensional stream function. 

00 - subscript denoting freestream conditions. 
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INTRODUCTION 


The primary function of an aircraft's horizontal tailplane is to balance 
the pitching moment which arises when the lift force developed by the wing is 
not applied through the airplane's center of gravity. Control of the airplane 
in the vertical plane is accomplished by varying this balancing moment so that 
the wing seeks a new equilibrium lift. The size of the horizontal tail thus 
determines the allowable c.g. travel as well as the lift coefficient which can 
be generated at any time. One does not gain a higher degree of control capa- 
bility, however, without accepting at the same time a performance decrement 
resulting from the higher drag of the larger surface. 

Since the advent of the airplane, designers have pondered the problem of 
positioning and sizing the horizontal tail to achieve both good aircraft per- 
formance and acceptable longitudinal stability and control. Their problem is 
compounded by the passage of the wing through the air ahead of tail. This 
produces (1) a downward inclination or downwash in the air approaching the 
tail, and, if the tail is directly in the wing wake, (2) a lower velocity at 
the tail due to the wing drag force. If the tail is sized assuming the onset 
flow to be freestream but is instead immersed in the wing wake, it will be 
unable to generate the expected trim moment. Thus, to maintain control capa- 
bility one must either locate the tail such that it never encounters such 
momentum deficiencies or else increase its size proportional to the kinetic 
energy deficiency actually encountered.. The latter course entails a perfor- 
mance penalty. Because of the non-linear response of the tailplane to onset 
flow angle, it is important to know the local flow direction everywhere in 
the tail region as well as those areas where momentum deficiencies exist. 

One would like to predict the wake location and the velocity deficit in the 
region of the horizontal tail to facilitate the process of positioning and 
sizing. 

As part of its program to provide technical support to the light aircraft 
industry and to aeronautical education in general, the National Aeronautics 
and Space Administration funded the present study at North Carolina State 
University. This work seeks to provide rigorous, easy-to-use methods (which 
now means computer programs) to predict the downwash and velocity profiles in 
the vicinity of the horizontal tailplane as well as detailed descriptions of 
the methods used, their historical antecedents, and some worked out examples. 
The study is thus analogous to previous studies (Ref. 1,2, 3, 4) which treat 
riding and handling qualities, aircraft performance, and wing and body lift 
and drag in a similar fashion. 



As will become evident from the literature review, the procedures for 
estimating downwash angle and wake position are better-defined than' those for 
estimating wake velocity profiles. The downwash calculation is usually 
approached by replacing the wing with a series of horseshoe vortices repre- 
senting both the wing lift and the vorticity shed by the wing. The downwash 
is a result of the interaction of the flow due to these horseshoe vortices 
with the freestream flow. It is most intense in a curved plane emanating from 
the trailing edge of the wing. The locus of this plane is usually obtained by 
integrating the local downwash angle from the wing trailing edge to some down- 
stream point of interest. More rigor may be added to this method by permitting 
the trailing vortices to deflect both spanwise and vertically so that their 
paths become streamlines. These streamlines then define the wake centerline 
position. This approach permits both upwash ahead of the wing and downwash 
behind to be calculated with a greater degree of accuracy than does the more 
conventional straight vortex method. The downwash computer program presented 
in this report employs the deflected vortex approach. 

Since the purpose of the present work is to provide a rigorous yet rela- 
tively simple method for determining the onset flow in the immediate neighbor- 
hood of the horizontal tail , it is not necessary to continue the calculation 
beyond this region to the point where the trailing vortices rollup into the 
f ami Mar vortical flow patterns observed far downstream of the aircraft. While 
the method employed herein can treat the inviscid flow aspects of this region 
sat isfactor I I y, other simpler methods are entirely adequate for developing the 
essential features of this "far field" problem; further, the method used here 
to determine the viscous dissipation present about the wake centerline is valid 
only in regions of quasi-pla nar flow, certainly not in vortical flow regions. 
For, this reason also, no attempt should be made to extend the calculation 
downstream beyond the region of the horizontal tail. 

The perceptive reader will note that no mention is made in the ensuing 
pages of the obvious contributions of the fuselage, both viscous and inviscid 
to the character of the flow field downstream of the wing. While the down- 
wash method used does include a treatment of the effects of fuselage inter- 
ference on the wing lift distribution, it does not include the effect of 
diminishing fuselage size as one moves aft on the downwash field nor the 
effects of the rapidly growing boundary layer on the aft fuselage. These 
were regarded as small perturbations whose determination would be excessively 
difficult given presently available analytical tools; for this reason they 
are not included. The user must therefore consider this fact in applying the 
results to new design situations. 

Another contribution which has been neglected in the computation is the 
modification in the wing’s downwash field due to the upwash field of the 
horizontal tail. Since the upwash angles are inherently smaller than the 
downwash angles and since the tail lift is generally less than i/5 the wing 
lift, this contribution has been neglected. 
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The method used herein is based on the lifting line representation of 
the wing; it is therefore applicable only to situations where the wing is of 
moderate- to-h igh aspect ratio, unswept, and relatively lightly loaded (the 
fuselage does not cover a very significant portion of the wing area). Situ- 
ations where these conditions are not satisfied will usually demand a much 
more complex lifting surface treatment. 

Another consequence of the lifting line representat ion is that the induced 
field in the immediate neighborhood of the wing is distorted by the concentra- 
tion of the wing’s bound vortices along a single line instead of being distrib- 
uted over the whole surface. Thus the downwash and upwash angles in this 
region may be excessive while at greater distances they will be indistinguish- 
able from those obtained from a lifting surface representation. Since, the 
region of interest in light aircraft is located more than 1.5 wing chords 
downstream of the wing’s trailing edge, a lifting line representat i on . shou I d 
be quite adequate. 

The variation in flow velocity as one moves vertically across the wing 
plane 'of symmetry is usually approximated by empirical relations derived from 
experimental data given the wing section drag coefficient and the distance 
behind the wing trailing edge. A more exact approach solves the two-dimen- 
sional boundary layer equations, using a finite difference technique, for the 
wake flow and then applies this solution along various streamlines displaced 
spanwise as given by the wake position procedure above. The result is a 
quas i -three-d i mens iona I solution. The starting conditions for the integration 
are the viscous velocity profiles on both upper and lower wing surfaces. This 
approach is used here to determine the vertical variation in flow magnitude. 

The work presented here begins with a two-part review of the pertinent 
literature dealing with the determination of (1) downwash angles and wake 
sheet location, and (2) velocity profiles in the wing wake. General theory 
and discussion of results sections are then presented for both the downwash 
(WASH) program and the velocity 'prof i le (WAKE) program. User instructions, 
program listings, and sample outputs of both programs are given in the 
append i ces. 
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LITERATURE REVIEW 




DOWNWASH AND WAKE POSITION 


In order to estimate the lift and drag forces experienced by an aircraft' 
horizontal tail, it is imperative that the magnitude and direction of the loca 
onset flow be known. The disturbance produced by the wing is the major con- 
tributor to non-uniform horizontal tail onset flows. Dwinnell (Ref. 5) points 
out that a wing generating a lifting force must have a pressure differential 
from the lower to the upper surface with the higher pressure on the lower sur- 
face. At the wing tips the air tends to flow laterally from the lower to the 
upper surface in an attempt to remove this differential. This lateral or 
spanwise flow is from wing tip to wPng root on top and from wing root to wing 
tip on the bottom (see total velocity vectors in Figure 1). The air immedi- 
ately aft of the wing tra i I i ng edge is given a swirling mot ion due to the 
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€ = downwash angle 


high pressure region 



upper surface 

lower surface 
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leading edge 
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trailing edge 


Figure i. Side and planform views of flow over a wing 

(bound and trailing vortex patterns not shown 
in sketch). 

shearing of the transverse velocity components on the upper and lower surfaces 
These swirls or vortices, more predominant at the wing tips, are referred to 
as trailing vortices, and they tend to give the air behind the wing an average 
downward inclination, usually called downwash (see Figure t for definition of 
downwash angle). The spanwise variation of the downwash is a function of the 
wing planform shape. 
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The circulation theory developed in the late 1800's and early 1900's to 
account for lift In an inviscid fluid proved to be very useful in estimating 
downwash. In 1878 Lord Rayleigh showed that the swerving flight of a "cut" 
tennis ball could be explained best in general terms by comparing it to the 
case of a cylinder placed in an inviscid uniform stream. By Imposing a cir- 
culatory flow upon the cylinder, the cylinder developed a force normal to the 
direction of the uniform stream (a lift force), directly proportional to the 
strength of the circulatory flow. In 1910 Joukowski showed that when a 
cylindrical body of arbitrary cross section moves with a velocity V in a fluid 
whose density is p and there is a circulation of magnitude r around it, a 
force normal to the velocity is produced equal to the product pVr per unit 
length of the cylinder. Joukowski made the assumption that the flow must 
leave the sharp trailing edge of airfoil-like bodies smoothly. Using this 
hypothesis, usually referred to as the Kutta-Joukowski condition because Kutta 
reached the same assumption independently somewhat earlier, the problem of 
lift became purely mathematical: one has only to determine the amount of 

circulation so that for zero vertex angle at the trailing edge, the velocity 
of the flow leaving the upper surface is equal to that leaving the lower sur- 
face. For finite vertex angles the trailing edge must be a stagnation point. 

Although the work of Rayleigh and Joukowski was concerned with wings 
having infinite aspect ratios (no spanwise flow), the circulatory flow idea 
was extended to wings with large aspect ratios and little or no sweep by 
Prandtl in his now well-known Lifting-Line Theory, von Karman (Ref. 6) 
summarizes Prandtl's theory in the following way: (a) the wing is replaced 

by a lifting line perpendicular to the flight direction; (b) the lifting line 
is assumed to consist of a bound vortex with circulation variable in order to 
account for the fact that the lift may change along the span; (c) in accor- 
dance with the change in circulation along the span, free vortices are born 
and extend downstream; however, (d) the flow produced by the vortex system is 
considered as a small perturbation of the fundamental stream relative to the 
wing, and therefore (e) it is assumed that the free vortices approximately 
follow the original direction of the streamlines parallel and opposite to 
the flight direction; (f) the flow in the immediate neighborhood of a wing 
section is determined by the two-dimensional solution given by Kutta and 
Joukowski . 

In the lifting line theory the wing surface is replaced by a line vortex 
of varying strength, called the bound vortex. A sheet of trailing vortices 
extends from the lifting line downstream to infinity: The perturbation in the 

flow field due to an infinitesimal segment of each vortex is found from the 
Biot-Savart Law (expressed here following Reference 7 in Vector notation). 


- r d£ x a 

w = — =— 

4tt | - | 3 
a 


( 1 ) 
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where w = induced velocity vector at a point in the flow field 

r_ - strength of the vortex whose segment is being considered 
da = vector having the direction and length of the infinitesimal 
vortex filament 

a = the vector from the vortex filament to the point 
the local value of r on the bound vortex is determined from 


Vc c £ ( y } 

r<y) = f 


( 2 ) 


where C £ is the section lift coefficient of the wing 
V is the free stream velocity 

c is the wing chord 

The strength of a particular trailing vortex is determined by the change in 
the r of the bound vortex which occurred as one moved outboard from the 
spanwise position where the previous trailing vortex was shed. All of the 
vorticity on the wing must eventually be shed downstream. The flow direc- 
tion in the field downstream (or upstream for that matter) of the lifting 
line is then a sum of the contributions of all the trailing vortices, the 
bound vortex, and the free stream. The perturbation in the free stream 
velocity field resulting from these contributions is most intense in a sheet 
which for positive lift moves downward and to the rear from the lifting line, 
in the physical situation which this model represents this sheet is essentially 
the center of the wing wake. From equation (2) it is evident that the greater 

the wing lift the greater the vortex strength; thus it is reasonable to expect 

that the greater the wing lift the greater will be the downward deflection of 
the wake centerline as it moves aft. 

Attention will now be directed to ways of evaluating the r distribution 
over the wing since it is essential for downwash calculations. Reference 8 
notes that the main problem of lifting line theory is the determination of the 
spanwise distribution of air loads on wings from the fundamental equation of 
the theory. Its simpliest form relates the absolute and effective angles of 
attack. 


a = a - a. (3) 

a o t 

where a = absolute angle of attack 
0 

a Q = effective angle of attack 
a. = induced angle of attack. 
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The effective angle of attack may be expressed in terms of r by writing the 
section lift equation as 


pvr = pv 2 


cm a /2 
o o 


or 


2r 

a o m Vc 
o 


where p = density of air 

m Q = section lift curve slope. 

Now the induced angle of attack may be written as - tan ^ (-^0 % - w/V and since 
w can be expressed in terms of r, then the only unknown in the fundamental 
equation is r. Using the Biot-Savart Law the downwash velocity w at the 

/q ' 

point y on the bound vortex, induced by a vortex behind the wing at the point 
y may be written as 


w 


y y 

7 o 7 


dr i 

4tt y - y 
7 o 7 


(4) 


Replacing dr by (dr/dy)dy 
all the trailing vortices 
b/2 (b is the wing span). 


and noting that the downwash at 
obtained by integrating w 


may be 
then 


V 


resulting from 
from - b/2 to 



+ 


b/2 



-b/2 


dr/dy 

y G - y 


dy . 


(5) 


The fundamental equation can now be written in the form 

b/2 


(a ) = (- 


2r 


V / x -) + 

ay m Vc y 4 ttV 

' r\ f r\ 


i- r 

ttV J 


-b/2 


dr/dy 

y„ - y 


dy 


( 6 ) 


where the subscript refers to the variable evaluated at the point y . 

Probably the best known procedure for solving the fundamental equation is 
given by Glauert in Reference 9. The first step is to invoke the following 
subst i tut ions : 


y = b/2 cos 0 


( 7 ) 
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and 


r = Vb y A sin n0 . (8) 

_ i n 
n=1 


The sine series rather than the cosine series is chosen for r because the 
boundary conditions of r = 0 at both wing tips (6 = 0 and 0 = tt) are auto- 
matically satisfied for ali values of A . The constants A are found by 
determining the values necessary to force the flow over the wing to be tangent 
to the mean camber surface. It can then be shown that the fundamental equation 
can be written as 


CO 

y CA sin n0(yn + sin 0)H = -5- a sin 0 (9) 

. n z a 

n=1 


where y = m c/(4b). 
o 

The above equation is valid for any spanwise station (0 between 0 and it), and 
in general y and a will be functions of 0. The problem has been reduced to 
one of determining the values of the coefficients A which will satisfy the 
above equation at every desired value of 0 along the span of the wing in 
question. If a 5 term series is desired for T then equation (9) must be 
applied at 5 spanwise stations yielding 5 equations In the 5 unknown A T s. 

■Once the A 's have been obtained then Equation 8 will yield r for different 
spanwise positions. 

As an alternative to solving Equation 9 at several spanwise stations, 
the computer programs given in Reference 4 are quite useful. Programs are 
presented to (1) calculate the lift, drag, and moment coefficients of conven- 
tional airfoils (two-dimensional) and (2) given the wing planform shape and 
airfoil character ist i cs for the root and tip sections, extend these two- 
dimensional characteristics to three dimensions. By using a lifting line 
theory, the three-dimensional extension yields local lift, drag, and moment 
coefficient information at several spanwise stations aiong the wing span. 

The section lift coefficients may be used directly in Equation 2 to compute 
the strength of the bound vortex. 

With the foregoing as background, attention will now turn to literature 
specifically treating the calculation of downwash. Two 1939 NACA Reports have 
become classics in the field: NACA TR-648 (Ref. 10) and NACA TR-651 (Ref. 7). 

The latter presents methods for predicting important wake parameters in terms 
of the distance behind the airfoil trailing edge and the profile drag coeffi- 
cient for plain and flapped wings. The methods were based on a correlation of 
extensive experimental measurements of downwash angles and wake characteristics. 
Downwash was computed in the wing chord plane using the Biot-Savart equation 
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operating on a vortex system made up of a bound vortex at the wing quarter- 
chord and a planar vortex sheet extending straight behind the quarter-chord 
line to infinity. It is pointed out that the wing wake coincides with the 
trailing vortex sheet because of their common origin at the trailing edge and 
their equal freedom to move in the induced velocity field behind the wing. 

The vertical displacement of the vortex sheet can therefore be approximated 
by integrating the down wash angle in the streamwise direction from the trail- 
ing edge to some specified distance downstream. Silverstein, Katzoff, and 
Bullivant used the idea of superposition to obtain maps of constant downwash 
angles behind the wing. First, using a planar vortex system, downwash angles 
were calculated at several points in the plane of symmetry to yield a contour 
map similar to that of Figure 2. This map was then distorted in the vertical 
plane to account for the vertical displacement of the vortex sheet. The 
distorted contour map is also shown in Figure 2. Using this procedure, which 
neglects spanwise curvature of the wake sheet as well as vortex roll-up effects, 
agreement with experimental data is good at distances typical of to the tail 
plane locations. 

In NACA TR-648 (Ref. 10) Silverstein and Katzoff prepared downwash charts 
to make the methods presented in TR-651 more readily useable in design 
activities. Included in the report are curves of the lift and drag of flapped 
airfoil sections and charts for finding the contribution of the different 
types of flaps to the total wing lift. The downwash charts are for elliptical 
wings and wings of taper ratios 0.1, 0.2, 0.3, and 0.5, with aspect ratios of 
6, 9 and 12, having flaps covering 0, 40, 70 and 100 percent of the wing span. 
These charts may be used without a prior knowledge of the material in TR-651. 
Although written some 35 years ago, these reports are still the basis of 
current downwash predictions for light aircraft. 

The interest in high speed flight led to some new downwash investigations 
in the late 40’s and early 50’s. In Reference 11 a study was made in the 
Langley 7- by 10-foot wind tunnel to determine downwash angles behind wings of 
various sweep angles and aspect ratios. Although the emphasis was on sweep 
effects, downwash characteri st i cs were also measured for some unswept wings. 

The results are displayed graphically as curves of constant downwash angles. 

In 1954 an experimental investigation was conducted by Tolhurst (Ref. 12) to 
define the downwash field and vortex sheet shape behind a large-scale wing- 
fuselage combination incorporating a 63° swept-back wing of aspect ratio 3.5. 
Data is presented for 9 vertica I -transverse planes behind the wing for 3 
angles of attack at a Reynolds number of 6.1 million. A comparison of these 
experimental results with theory showed that an acceptable approximation to 
the downwash distribution within the distance surveyed could be obtained If 
the vortex sheet were assumed to be flat and the spanwise vorticity distri- 
bution were assumed to be that of the wing. Rogers (Ref. 13) studied the 
application of two-dimensional vortex theory to the prediction of flow fields 
behind wings of wing-body combinations at subsonic and supersonic speeds. The 
wing trailing vortex was represented initially by line vortices distributed to 
approximate the spanwise distribution of circulation along the trailing edge 
of the exposed wing panels; however, a bound wing vortex was not included. 
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Figure 2. Example of downwash calculation procedure (Ref. 7) in the 
plane of symmetry for a USA 415 airfoil at = 1.175. 
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The afterbody was represented by correspond i ng Image vortices within the body. 
Two-dimensional line-vortex theory was then used to compute the induced 
velocities at each vortex, and the resulting displacement of each vortex was 
calculated by means of a numerical stepwise integration procedure. The dis- 
placement of each vortex independently permits inclusion of both the roll-up 
phenomenon and the spanwise wake-sheet curvature. Rogers found that this 
method predicted vortex paths at high speeds with reasonable facility and 
accuracy. 

In 1970 Butter and Hancock presented a method (Ref. 14) for calculating 
the shape of the trailing vortex system behind a swept wing at low speed. 

The authors were concerned about the two-dimensional approaches which omit the 
main wing bound vortex; they believed that when C. is large the wing circula- 
tion effects are significant. A method was therefore proposed which would 
include the effect of the main wing bound vorticity and the three-dimensional 
pattern of trailing vorticity. The wing is replaced by a single lifting-line 
bound vortex which represents the circulation distribution around the wing 
while the trailing sheet is approximated by a number of discrete line vortices 
Starting with a nearly planar system (_[_.e_. , bound vorticity on the wing and 
trailing vorticity streamwise as in Figure 3) and using a step-by-step process 
working downstream aft of the wing trailing edge, the vertical and spanwise 
displacements of the trailing sheet are determined. To calculate the position 
of the ith vortex at station 1, fix behind the 3/4-chord line (station 0), the 
sidewash and downwash velocities, (v. and w.) due to all the other vortices 
are calculated at 6x/ 2. The change in vertical and spanwise position (Sz and 
6y) for the' ith vortex from station 0 to station 1 is then given by 


Sy = v.5x/V and <5z = w ,5'x/V 

7 i i 


( 10 ) 


This procedure is used to obtain the new position of each vortex at station 1. 
Downstream of station 1 each vortex extends to infinity using the y and z 
coordinates at station 1. Once the locations are established at station 1, 
the remainder of the stations downstream may be analyzed in a similar manner. 
Butter and Hancock found that the principal factors governing the accuracy of 
the results are the downstream step-size and the number and position of the 
vortices chosen to represent the trailing vortex sheet. They note that If too 
many trailing vortices are taken numerical difficulties may arise; induced 
velocities calculated at a point between two such vortices will be given by 
the difference of two large velocities of opposite sign. When compared with 
experimental results, the downwash angles were found to be too large near the 
trailing edge due to the concentrated nature of the bound vortex; however, 
further downstream agreement was good. 

In Reference 15 Labrujere presents a method for calculating the vortex 
location behind a thick lifting wing. The flow around a wing in uniform onset 
flow is determined with the aid of a source distribution on the surface of the 
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+ for the method presented in Referee 

Figure 3. 




wing and a vortex distribution on its camber surface and wake. The wing and 
camber surface are divided into a number of strips whose edges are formed by 
lines y = constant and which are each subdivided in a number of quadrilateral 
panels. Along the periphery of each camber panel a line vortex is located 
(see Figure 4). The strengths of these "singularities” are determined by 
imposing the boundary condition of tangential flow to the wing surface at a 
number of surface points and the condition of tangential flow to the trailing 

vortex sheet at a number of "Kutta points" noted in Figure 4. Since the 

position of the vortex sheet and the strength of the singularities are related 
"non-l inear ly", an iterative procedure is required in order to determine the 
shape of the sheet as well as the singularity distributions for a given wing 
and at a given angle of attack. Although only a limited comparison with 
experimental results was made using this method, the agreement was good. 

Comparing Labrujere’s technique with that of the previous reference, 
there are three reasons why the Butter and Hancock method is preferred as the 
basis for further work: (1) Labrujere combines the problems of the flow over 
the wing and the flow in the wake, solving them simultaneously; for a large 

number of panels both on the wing and in the wake this requires both more 

computer storage and a longer execution time than is required for the two 
problems independently (necessary in the Butter and Hancock approach). (2) 

The iterative nature of Labrujere's method may necessitate several executions 
unless the initial estimate of the wake shape Is accurate. (3) Although Butter 
and Hancock's method may be used iteratively to improve accuracy, they found 
that for step sizes of a quarter semispan and smaller, subsequent iteration 
did not significantly affect the position of the trailing vortices. 

In the last few years the interest in trailing vortex sheets has turned 
to the far field problem. At large distances behind the wing, vortex roll-up 
has occurred, and the hazard of one aircraft flying into the vortex core of 
another aircraft has spurred interest in trying to analytically predict the 
shape and velocities of a rolled vortex pair as well as the dissipation of 
this vortex pair with time. References 16, 17, and 18 cover the topic very 
well and discuss the ■ pert i nent work which has been done to the present time. 
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WAKE VELOCITY PROFILE 


While an inviscid theory Is acceptable for predicting the location and 
inclination of the wing wake, the estimation of actual velocities requires a 
viscous analysis. The velocity deficit behind a wing is a function of the 
wing drag, a viscous phenomenon caused by the shear forces which the wing 
introduces in the fluid. The only really correct way to attack the problem 
is to solve the full, time-dependent Navier-Stokes equations of motion in 
three-dimensions over a substantial region. With present computers and com- 
putational techniques, computers will have to grow in storage size and speed 
by one or two orders of magnitude in order to obtain an accurate solution to 
the problem in the space of a few hours. 

Boundary layer theory, which has been a valuable tool in solving for the 
viscous flow over airfoils and even bodies at zero or small angles of attack 
(Refs. 19 and 4), has also been applied to the flow in jets and wakes. 

Boundary layer theory consists of simplifying the Navier-Stokes equations so 
that they are applicable in the viscous flow near the surface of a body or in 
a wake. For the flow over a body the region where viscous effects are dominant 
is usually very thin compared to the body size. The equations are simplified 
by estimating the order of magnitude of each term and neglecting those which 
are smal I compared with the other terms. As might be expected, these simplified 
equations are much easier to solve, yet they retain the correct physical 
description of the flow in most cases. The conventional way of approaching 
the wake problem is to assume that the flow downstream of the trailing edge 
can be described by the boundary layer equations which are then solved by a 
momentum integral or finite-difference method. The remainder of this section 
will review some of the pertinent literature on the determination of wake 
velocity p rof i I es . 

Two-dimensional wakes were probably first investigated by H. Schlichting 
in his thesis presented to Goettingen University in 1930. In his general text. 
Boundary Layer Theory (Ref. 20), he notes that a wake is formed behind a solid 
body which is being dragged through a fluid at rest or behind a solid body which 
has been placed in a stream of fluid. The velocities in a wake are smaller 
than those in the main stream and the velocity losses in the wake amount to a 
momentum loss in the fluid due the body drag. The spread of the wake increases 
as the distance from the body is increased, and the differences between the 
velocity in the wake and that outside become smaller. A wake is a free turbu- 
lent flow (not confined by solid walls) and therefore of boundary layer nature. 
This means that (1) the region of space in which a solution is being sought 
does not extend far in the transverse direction, as compared with the main 
flow direction, and (2) the transverse gradients are large. The problem can 
therefore be studied using the two-dimensional, steady, boundary layer 
equations given below: 

u + v = — —■ (x-momentum equation) (11) 

3x 3y p 3y M 
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au 9 v 

9x 9y 


0 


(continuity equation) 


( 12 ) 


where x and u are the coordinate and velocity along the wake centerline, y and 
v are the coordinate and velocity normal to the wake centerline, p is the 
density of the fluid, and t is the turbulent shearing stress (see Figure 5). 


velocity profiles 



Figure 5. Sample velocity profiles in an airfoil wake. 


The pressure gradient term has been dropped from the x-momentum equation because, 
from a certain distance behind the body on downstream, the pressure remains 
essentially constant in the y-direction. This assumption is not necessarily 
true near the airfoil trailing edge. To integrate Equations 11 the turbulent 
shearing stress must be expressed in terms of the main flow parameters. One 
of the simpliest and most widely used wake flow models for t is that given by 
Prandtl in Equation 13. 


T 



pcb ( u 


max 


u . 
mm 



(13) 


where b denotes the width of the wake, c is an empirical constant, and e is 
the virtual kinematic viscosity or eddy viscosity, assumed constant over the 
whole wake width. Schlichting points out that experience has shown that for 
two-d imens iona I turbulent wakes, the wake width increases as /x while the 
centerline velocity decreases as ]//x. 

As noted in Reference 20 Goertler and Reichardt investigated the two- 
dimensional wake behind a circular cylinder. Goertler found that if he 
assumed : 
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(1) u 1 = U - u which at large distances from the body is small compared 

to U, the freest ream velocity, 

(2) the drag on the cylinder is obtained by applying the momentum theorem 

far downstream so that 


y=+oo 


D = hp 


I 


u 1 ( U - 


u 1 )dy 


(14) 


where D is the drag on the cylinder of height or length h, 

(3) that the second term of Equation 11 is small compared to the first, 

then, by integrating the simplified version of Equation 11 


u 4/a 


uc 0 d 


1/2 


( ^- ,/2 <-n 2 /4) 

PH ® 

C D d 


where : 


= cylinder drag coefficient 
e Q = constant eddy viscosity 
d = cylinder diameter 

x = distance downstream behind the cylinder 
n = nond imens iona l i zed y-coordinate = y 



(15) 


Upon comparison with experimental measurements made by Schlichting, Goertler 
evaluated the constant in Equation 13 with the relation 


e 

o 


(0.047)(2b 1/2 u ] ) 
max 


where: b i/2 = ha ^ ^ he wa ^ e width halfway between and Uj at some 

x-station max min 

u 1 = the maximum value of u^ in the wake at that x-station. 
max 


Using similar assumptions and a slightly different eddy viscosity model, 
Reichardt arrived at the following expression for u^ : 


/To 

18(. 18) 


(-*-> 

C D d 


^ ^ / 3/2\^ 

\1 - (y/b)^/ 


(16) 
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Schiichting (Ref. 20) notes that the two solutions for the wake velocity 
actually do not differ much one from the other. 

While the above work was concerned with the wake of a circular cylinder, 
in 1939 Silverstein, Katzoff, and Bullivant (Ref. 7) generalized the available 
experimental results so that velocity profile and wake width information 
behind airfoils could be predicted in terms of the downstream distance and the 
profile drag coefficient. The profile drag of an airfoil may be approximately 
equated to the loss of momentum in the wake by the relation: 


d 


o 



u) dy 


(17) 


where d c is the drag per unit span. At distances behind the wing comparable 
with tail-plane location, where the static pressure in the wake has reached 
that of freestream, experimental investigations have shown that Equation 17 
may be approximated in coefficient form by the more elementary expression: 




( 1 - 3-) dy 
q 

^o 


(18) 


where. 


c d = section profile drag coefficient 

q/q° = ratio of dynamic pressure in the wake to that in the free- 
^ o , 

stream 

b = wake w i dth 
c = wing chord. 


The wake may be completely described by the width b, the loss of dynamic 
pressure at the wake center Aq, and the shape of the wake profi le. As an 
aid in generalizing the results, the following nond imens i ona I ratios were 
adopted : 

n = M = (| _ q/q Q ) 

^O 

_L , wake width . __ b__ 

^ 2 wing chord 2c 

distance behind trailing edge of wing _ l_ 

^ wing chord c 
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Using Prandtl f s proportiona I i+y relations 


1/2 


5 * C 1/2 c. 


and 


n 'v 


. 1/2 


(19) 


the available experimental data was used to find the proport tona I i ty constants 
which should be used in Equations 19. It was found that: 


C 


0.68 



(5 + 0.15) 


1/2 


( 20 ) 


t > 


2.42 



o 

£ + 0.3 


( 21 ) 


fit the data very well. The equation for n was found to vary with 1/£ rather 
than 1 /S 1 ^ as mentioned previously. The authors point out that discrepancy 
is probably associated with the high values of n near the trailing edge of the 
experimental data used for correlation; in Prandtl’s formulation n was assumed 
small. For a complete wake descri pt ion, the shape of the profile must also be 
defined. The profile shape was found to be closely approximated by either of 
the following empirical relations: 


n' 


n 


D 



1.75 2 


( 22 ) 


or 


ai = 

n 


2 , TTg * 
cos (^) 


(23) 


where n f = relative loss in dynamic pressure 

= distance from wake centerline (?’ < O . 

The profiles given above are assumed symmetric about £’ = 0 and for E, values 
in the vicinity of the tail plane this is not a bad approximation even for 
unsymmetric airfoils. 

In TR-648 (Ref. 10) Silverstein and Katzoff made the methods presented in 
TR-651 easily adaptable to design by preparing charts for determining the 
velocity profile at the horizontal tail. Over the years both TR’s have been 
very valuable in predicting tail plane velocity profiles for light aircraft. 
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In 1947 an investigation of free turbulent mixing was presented by 
Liepmann and Laufer in Reference 21. The investigation examined the various 
ways turbulence had been modeled up to that time and discussed the physical 
significance of the models using experimental measurements of the mean 
velocity, the intensity and scale of turbulent fluctuations, and the tur- 
bulent shear in a two-dimensional mixing zone. The authors decided that the 
assumptions of constant mixing length and constant exchange coefficients (as 
in Equation 13) across the wake width were, strictly speaking, incorrect. 
Measurements of the field fluctuating velocities in a turbulent mixing zone 
show that both the mixing length and exchange coefficient vary across the 
mixing region, contrary to the assumptions made by Prandtl and others. They 
concluded that if only the mean-velocity distribution of a turbulent flow 
problem is desired, it usually can be obtained by dimensional considerations, 
assuming a reasonable curve (_e.il. > an error function in the case of a mixing 
region), and determining one constant from experiment. (This agrees with the 
methods already presented from References 20 and 7). Also, for an understand- 
ing of turbulent flow in general, a study of various mean-velocity distribu- 
tions is of little use; the field, of turbulent fluctuations must be studied 
in detail probably requiring a knowledge of the diffusion of turbulent energy. 
The constant exchange coefficient as presented in Equation 13 can therefore 
be used in a turbulence model with good success if only mean velocity profiles 
are desired; but for any deeper understanding of turbulence, the micro- 
structure must be investigated more fully. 


The subject of eddy viscosity models is treated to some extent in 
References 22 and 23. Elassar and Pardolfini, using a multilevel linear 
difference scheme examined four eddy viscosity models and evaluated the con- 
stants in these models by comparing calculated velocity profiles with experi- 
mental data in the similarity region. The thrust of the report is directed 
toward high-speed compressible flow, and all the models are more complex than 
Equation 13. In Reference 23 similarity solutions of the boundary layer 
equations in wakes with tailored pressure gradients were obtained. Experi- 
mental data was used to obtain a value for the constant eddy viscosity and the 
agreement between theory and experiment is described as good. While both of 
these reports are valuable in a general nature, it would be desirable to have 
viscosity constants evaluated using data obtained in the wake behind an air- 
foil. References 24 and 25 indicate that boundary layer solution procedures 
are being applied to high-speed turbulent mixing layers using more elaborate 
viscosity models than that given in Equation 13. These same techniques should 
be readily adaptable to low speed turbulent mixing and could be very useful 
If computer storage requirements and execution times are not excessive. 


in order to improve upon the turbulent wake prediction methods, more exper 
imental data are necessary for wake profiles at several stations behind an 
airfoil rather than at the tail-plane location only. References 26 and 27 are 
very helpful in this respect in that they present two-dimensional velocity 
profiles in the wakes of symmetrical Joukowski and Piercy airfoils respectively 
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An 11.8 percent thick Joukowskl airfoil was tested at angles of attack of 0 
and 6 degrees for a Reynolds number of 0.42 million. Velocity. prof i les were 
measured at stations of 0.0005, 0.0025, 0.01, 0.25, and 0.5 chord lengths 
behind the airfoil trailing edge. A similar investigation on the 12 percent 
thick Piercy airfoil, with wires on the leading edge to force transition, 
yielded velocity profiles at 0.0005, 0.0025, 0.01, 0.25, and 0.4 chord lengths. 
In addition to the profiles, the drag coefficient and other boundary layer 
information was also measured. 

A recent report by Goradia and Li I ley (Ref. 28) also present some experi- 
mental data for profile measurements behind a modified NACA 0015 airfoil at 
several angles of attack. The report describes a procedure developed by the 
authors to improve two-dimensional drag predictions for single and multi- 
element airfoils by describing the airfoil wake using regions of similarity. 
Although correlations are given for only two or three airfoils and the computer 
program is not provided, the procedure seems to have improved the drag predic- 
tions given in CR-1843 (Ref. 29). Experimental data, taken in connection with 
this work is presented in both tabular and graphical form; however, the 
tabular data, particularly at stations near the airfoil trailing edge, 
indicate profiles which fail to return to freestream velocity at the upper and 
lower edge. 

While the data given in the three references cited above are valuable, 
they were taken in the wake of symmetrical airfoils. More data is obviously 
needed, part i cu I ari 1 y for unsymmetric airfoils, but the lack of data is due 
in part to the high cost of wind tunnel testing and the difficulty of measur- 
ing wake velocities accurately. 

A review of the literature revealed several papers concerning the effect 
of flow curvature on the formulation and solution of the two-dimensional 
boundary layer equations. Murphy (Ref. 30) formulated the boundary layer 
equations with curvature and indicated that even for laminar flow these 
equations are at their simplest more complex than the no-curvature case. 
Williams (Ref. 31) obtained a solution for a free jet flowing off a step only 
by assuming a curvature function which made the boundary layer equations 
similar (the streamwise coordinate did not appear explicitly). One encounters 
similar increases in solution complexity when treating three-dimensional 
boundary layers. Although the problem can be formulated, the computer storage 
requirements and execution times required with present machines and methods 
are excessive. Add to this the complexity of curvature and it is understand- 
able why attention is still directed toward solving the simplified problem. 
Solution to such a problem, hopefully, preserved the qualitative aspects of 
the physical situation even though it cannot duplicate it quantitatively in 
a I I respects. 
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PROGRAM FOR DOWNWASH AND 
WAKE POSITION CALCULATIONS 
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GENERAL THEORY 


For many years the vertical displacement of the wake centerline location 
and downwash behind wings has been calculated using a bound vortex and a 
finite number of trailing vortices. The spanwise curvature of the trailing 
vortices and their roll-up may be included in the calculation by allowing the 
vortices themselves to be deflected in both the vertical and spanwise direc- 
tions (Ref. 13 and 14). Analyses have also been conducted by representing the 
wing and vortex sheet as a single vortex lattice and using an iterative 
solution procedure with boundary conditions of no flow normal to the vortex 
panels to determine the wake location (Ref. 15). Because of the paneled wake 
procedure (lattice) achieves only a slight improvement in wake location at a 
large increase in computer cost , the WASH program presented here uses a 
method similar to those given in References 13 and 14. For simplification 
the x, y, and z coordinates are non-d imens iona I i zed by the wing semispan, and 
a unity freestream velocity is assumed since only ratios of- induced velocities 
to freestream velocity are required to estimate the downwash and sidewash 
ang I es. 

The first step in computing the downwash characteristics and wake 
position is that of representing the wing by a bound vortex at the wing 
quartei — chord and by a trailing vortex system running aft from the quarter- 
chord line. The symmetric ha If -wing is placed in the coordinate system shown 
in Figure 6. The root chord I i ne is in the xz-p I ane at y = 0, with the origin 
at the trailing of the wing root. The root chordline has an incidence, a, 
with the negative x-axis. 



Figure 6. Definition of coordinate system and wing placement used in 
WASH program. 
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The analysis will be limited to unswept wings; therefore, the x and z 
coordinates of the wing quarter-chord are given by 


3c cos a 
x b 2b 


and 


3c sin ot 
Z b 2b 


where ^ x b ,,y, b ,z b^ = "*^ e coorc * i nates of the bound vortex normalized by the 

semi span 
b = wing span 
c = wing root chord 
a = wing angle of attack. 

The strength of the bound vortex is calculated using Equation 2 once the 
section lift coefficient distribution is known. Although the WASH program 
was designed to use the output of the three dimensional aerodynamic charac- 
teristics program given in Reference 4, any procedure may be used to calculate 
the section lift coefficients. The lift coefficients at 10 different 
normalized spanwise stations must be read into the WASH program. The speci- 
fication for both is ordered from positive wing tip to wing root (see Figure 
7). 


wing 

root x denotes specified spanwise stations 

L * K X X X X X x- 

10 98 765432 


wing 

tip 


I 


* 


I 


y 


Figure 7. Ordering scheme for spanwise station and section lift 
spec i f i cat i on . 


Recall of course that the wing and its flow field is being represented by a 
system of horseshoe vortices. The initial streamwise position of the various 
trailing vortices is defined by assuming the vortices are shed at stations 
midway between the 10 specified y values on each semi-span, j_.e^. , on each half 
of the lifting line. In carrying out the computations, advantage can be taken 
of the assumed symmetry in the loading and the geometry. In this fashion, the 
section lift coefficients as well as the input y stations and the positions of 
the shed vortices for the negative y axis can be defined. A schematic of the 


* The root value for y cannot be zero; it may be a very sma I I number, but 
because of the technique used for estimating the bound vortex contribution, 
the value must be non-zero. 
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20 stations at which the section lift coefficients are specified and the 18 
shed vortex stations for the entire wing is given in Figure 8. 


x denotes specified spanwise stations for complete wing 
I denotes shed vortex stations for complete wing 


wing 

tip 

| I 2 3 


wing 

root 


. . . 16 17 


wing 

tip 

18 


I 2 3 . . . 


. . . 18 19 20 


Figure 8. Reordered specified spanwise stations and shed vortex 
stations for complete wing. 


Once the bound vortex strength is calculated using Equation 1 at the specified 
y stations, the strength of the shed vortices is computed by the change in r 
from one station to the next. Thus, for the right half of the wing the trail- 
ing vortex system is similar to that given in Figure 9. 



Figure 9. Example of indexing procedure used for the trailing vortex system 
on the right half of the wing with a typical r(y) distribution. 
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Each initial trailing vortex is made up of two straight line segments, one 
from the bound vortex to the trailing edge and another from the trailing edge 
to infinity in the x-di recti on parallel to the freest ream (see Figure 10). 



Figure 10. Example of initial approximation of trailing vortex system. 


The net velocity induced by the vortex system at a point in the vicinity 
of the wing is computed by summing the contributions from the bound vortex and 
every straight line segment of each trailing vortex. The velocity induced at 
a point due to a vortex is defined by the Biot Savart law (Equation 2). An 
alternate form of this equation for straight line filaments is given in 
Reference 5 as 


r sin 9 dj, 

4t r 2 

r 


(24) 


where 


dW = total velocity induced at a particular point by a filament d£ 
r = strength of the vortex filament 

r = distance to increment of filament at the point the velocity is 
induced 

0 = angle between the length r and the filament. 
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The contribution of the bound vortex is calculated by integrating Equation 24 
from the negative wing tip to the positive tip. Expressing sin 0 and r in 
terms of the x, y, and z coordinates, and noting that d£ becomes dy, the 
induced velocity at a point (x, y, z) due to the bound vortex is written as 


w b ■ 


i 

; 


r(y o ) b/2 \ 


(x - x.) 2 + (z - z.) 2 dy 
b b 7 o 


Ii 4tt[(x - x ) + (y - y ) + (z - z.) 1 


2-i3/2 


(25) 


where all coordinates are non-d imens iona I i zed by the semi span. Equation 25 
may be solved numerically, or if a functional form is known' for r(y Q ), it may 
be integrated in closed form. In order to use the closed form technique, 
two successive specified values of r(y) for the bound vortex were fitted using 
the cubic spline given in Reference 33 matching first and second derivatives 
at the end points. The cubic spline fit is of the form 

rty) = A,/ + A 2j y 2 + A 3j y + A 4j (26) 


where y. < y < y. + ] , j = 1, 2,..., 19. 

Note that for each of the 19 intervals the cubic coefficients depend on the 
interval (yj, y- + ^) in which the spanwise station lies, but they are constant 
in a given interval. It is this variation of the coefficients that gives the 
spline fit its remarkable curve fitting properties. Using the cubic represen- 
tation Equation 25 is integrated (Ref. 34) in each of the 19 intervals yield- 
ing the total velocity induced by the bound vortex at the point (x, y, z). 

This total velocity is divided into a freestream component, u, and a vertical 
or downwash component, w, by the relation given below: 


u b = v* - x b )/ V (x • v 2 + <z - z b )Z 
w b = W b (z ' z b ,/ V (x ■ v 2 + <z * z b )2 


(27) 


w and W are defined as positive in the negative z direction, and u is positive 
in the positive x direction. Since the bound vortex is always parallel to the 
y-axis, there is no si dewash component in W^. The streamwise component u^ is 
usually very small compared with the freestream velocity. 
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The velocity induced by a trailing vortex' is calculated by summing the 
effect of each straight line segment of the vortex. The total velocity 
induced at a point by a straight line vortex segment is also calculated using 
Equation 24; the evaluation is however, much simplier than the bound con- 
tribution because r along all segments of a trailing vortex is constant. 

Thus, expressing d£ and r in terms of 0, one has 



where W t = is the total velocity induced at a point (x, y, z) by a 

vortex of length AB, 

6^ and 0^ = angles between the length r and the filament at points 

A and B respectively, and h is the perpend i cu I ar distance 
from the point (see Figure 11). 



Figure 11. Nomenclature diagram for velocity induced by a straight-line 
vortex filament (Ref. 5). 


Carrying out the integration indicated in Equation 28 yields 


W 


t 


r 

4Trh 


( cos 0^ 


- cos 0 ) 


(29) 
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This total induced velocity must also be divided into its component parts: 
freestream, sidewash, and downwash velocities. As seen in Figure 12, the 
total induced velocity at P Q due to a line vortex from to - P 2 of strength r 
may be rewritten as 


T(cos 0^ + cos e 2 ) 
4ir b/2 d Q j sin 0^ 


(30) 



Figure 12. 


Schematic for velocity 


P to P 


2 * 


induced at P^ due to a vortex from 


The components are derived using the following relations to obtain the 
direction cosines: 

d 21 = V <X 2 ~ X 1 )2 + (y 2 " Y 1 )2 + (Z 2 " z 1 )2 

d oi = V (x o " x i )2 + (y o " y i )2 + (z o " z i )2 

d 02 = V (X 0 ' X 2 >2 + (y 0 " y 2 )2 + (Z 0 “ z 2 )2 


(31) 
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2, 2 * 

«o, 2 - 

d 02 >/<2d 21 d 01 

>j 

II 

CM 

<D 

cos 

^(d 

21 2 + 

d 02 2 - 

d 01 2)/(2d 21 d 02 

>] 

X 3 = 

x i + 

d 01 

cos 0 

1 <X 2 

- x,)/d 21 


y 3 = 

n + 

d 01 

cos 0 

1 (y 2 

- y,)/d 21 


N 

OJ 

II 

z i + 

d 01 

cos 0 

1 (Z 2 

- z,)/d 2I 


A = 

<y 2 ' 

- Vi 

><z 3 - 

v - 

(y 3 - v 1 )(z 2 “ 

z 1 } 


B = - C(x 2 - - z - (x^ - x^(z 2 - z^)D 

C = C(x 2 - x 1 )(y 3 - y ^ ) - (x 3 --x^)(y 2 - y^H . 

The components can therefore be written as 

u t = W t A//R 
v = W B//R 
w = W C//R 


(31 ) 
Cont. 


(32) 


where R = A 2 + B 2 + C 2 . 

In summary, the induced velocity components at a point are calculated 
due to U) the bound vortex using Equations 25 and 27 and (2) a straight 
line segment of a trailing vortex using Equations 30 and 32. The velocity 
due to the entire vortex system is computed by summing the contributions from 
the bound vortex and each segment of the 18 trailing vortices. 

The spanwise curvature in the trailing vortex system is modeled by repre- 
senting each vortex by a number of straight-line segments which are free to 
follow the local direction of the flow. The position of these vortex segments 
is determined using a stepwise downstream marching procedure which will now 
be discussed. 
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As previously mentioned, the first approximation of the vortex sheet 
consists of the !8 trailing vortices with 2 segments each, one segment from 
the bound vortex to the y-axis and another from the y-axis to x = 00 paral lel 
to the x-axis. Since the actual vortex system behind a wing is displaced in 
both the y and z directions, the following procedure is used to compute the 
position of each trailing vortex at a length Ax behind the trailing edge. 
Noting Figure 13 and considering vortex #10 at station #2, the v and w. 
induced velocity components due to all the other vortices are calculated 
halfway between stations #2 and #3 (x = Ax/2) using the y and z coordinates 
at station #2. The changes in spanwise and vertical position of vortex #10 
in moving from station #2 to #3 is approximated by Equations 10 as 


Ay = vAx 
Az = wAx . 



Figure 13. Schematic of downstream stepwise technique used in WASH 
program. 
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The position of the vortex at station # 2 was (0, y( vortex #10), 0) and the 
position at #3 becomes (Ax, y( vortex #10) + Ay, Az). A similar procedure is 
used to calculate the position of the other 17 trailing vortices at station 
#3. Downstream of station #3 the trailing vortices are assumed to go to 
infinity in the streamwise direction. After the position is determined at 
station #3 the procedure is reapplied at stations #4, #5, etc. until the 
vortex sheet is traced as far downstream as desired. The WASH program prints 
the coordinates of each vortex as well as the total vortex length from the 
trailing edge at each downstream station. The procedure discussed above is 
equivalent to that used in Reference 14 where it was found that for stepsizes 
as large as .25 semispans an accurate wake location was predicted. Butter 
and Hancock did note that numerical difficulties were sometimes encountered 
when more than 5 spanwise stations for the half-wing were used; however, the 
WASH program using double precision arithmetic utilizes 9 spanwise stations 
without numerical difficulty as long as the specified y stations are approx- 
imately evenly spaced. An example of a wake sheet calculated using the above 
procedure is given in Figure 14. The two views of this wake sheet were 
plotted using the PLOT program of Reference 4. The WASH program will punch a 
data set which may be read directly into the PLOT program; the data set 
requires only the addition of view cards at the rear defining the type of view 
desi red. 

Once the location of the trailing vortex sheet is determined then both 
downwash and upwash angle may be computed behind and ahead of the wing 
respectively using the bound vortex and the new vortex sheet. The WASH 
program is designed to calculate, at the user’s option, three basic sets of 
i.n format! on : 

(1) downwash and si dewash angles in a plane behind the wing at a con- 
stant x-station (a negative value for the x-station yields upwash 
angles ahead of the wing) indicating the variation with y and z, 

(2) downwash angles in the y = 0 plane of symmetry for 1.75 semi spans 
downstream with z varying from .5 semi spans above and below a user 
specified z-station, 

(3) upwash angles in the y = 0 plane of symmetry for 0.5 semi spans 
upstream with z varying from .25 semispans above and below the 
leading edge z coordinate. 

For each option presented above, a grid of coordinates, at which the 
downwash information is calculated, must be established. The f i rst option 
was originally designed to compute the downwash in any x = constant plane at 
all the grid points determined by y varying from 0 to 1.1 semispans in steps 
of .05 semispans with z ranging from .25 semispans above and below a specified 
z station in steps of .05 semispans. This required calculating downwash and 
sidewash angles at 253 points. Numerical difficulties were usually encountered 
at one or more of the grid points because of the grid points being located 
very near a, vortex segment. The vortex velocities by definition are unbounded 
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at the vortex center; thus, unrealistic velocities may be predicted if the 
grid point is too near a vortex segment. To correct this difficulty an 
attempt was made to determine an optimum distance to a filament such that the 
filament's contributions would be included or excluded depending on whether the 
distance to the filament was greater or less than the optimum. While this 
procedure was . successf u I in producing reasonable values of downwash angle, it 
does not guarantee a smooth variation in downwash In either the y or z direc- 
tion since the downwash at one of two closely-spaced points may include the 
contribution from a filament while the downwash computed at the' other point 
may not. 


To insure a smooth variation, the following procedure is used. The user- 
specified x-station at which the downwash information is required must be a 
multiple of the stepsize; for upwash calculations any negative x-station ahead 
of the wing leading edge may be used. The y and z coordinates for the 6 most 
inboard vortices on the positive y-axis are used to establish a new grid 
system. Using the y coordinate at one of the six vortices as constant, the 
total flow angularity (vector sum of downwash and sidewash angles) and the 
si dewash and downwash angles are calculated for z ranging from .25 semi spans 
above and below the z coordinate of that vortex. This procedure is also 
carried out for the other five most inboard vortices yielding a grid system 
similar to that given in Figure 15.. 



Figure 15. Example spacing for downwash calculations using NX0RY=1 
option (spanwise variation). 
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The numerical difficulties encountered because of- the contribution of the* 
vortex segment to the field in its immediate vicinity are eliminated by (1) 
excluding, for example, the contribution of vortex 11 when calculating the 
flow angles above and below the z coordinate of vortex 11, and (2) calculating 
downwash above only the six most inboard vortices which, for x-distances com- 
parable to horizontal tail locations, have not rolled-up and therefore do not 
pass above or below another vortex trace. The exclusion of vortex 11 when 
determining the flow angles above and below it will not affect the downwash; 
however, the si dewash will be in error. Since the downwash and not si dewash 
is the quantity which the program is designed to predict, the error in side- 
wash is tolerated. The option discussed above is produced when the input 
variable NXORY is specified as 1 and STAT is read as the desired x-station in 
semispans. A detailed description of these and the other necessary input 
variables is given in the User Instructions section of Appendix A. 

Option number two in the WASH program computes the dowriwash angles in the 
plane of symmetry, y = 0. The grid is specified by varying x from .25 to 1.75 
semispans behind the trailing edge in steps of .05, and z ranging from .5 
semispans above and below a midpoint z-va I ue specified as the variable STAT. 

As above, numerical difficulties were encountered in some cases using this 
option because the most inboard shed or trailing vortex on a half-wing may be 
too close to some of the grid points. To remedy the problem In general, the 
downwash angles were calculated using only the 8 most-outboard vortices on 
both sides of the y = 0 plane, and the WASH program is written in this fashion. 
However, if the user desires to include the effect of the inboard vortex, the 
DO LOOP parameter may be changed from 8 to 9 on card PLA 110 of the WASH 
program. Inasmuch as the computations are made in the plane of symmetry, only 
half the wing need be analyzed to compute the downwash, and the si dewash is 
identically zero. This downwash information is computed in the plane of 
symmetry when the variable NXORY = 2 and STAT is the desired midpoint value 
of z in semi spans. 

If NXORY = 3, then upwash information is calculated in the plane of 
symmetry. The grid points at which upwash is computed are established by x 
ranging from the wing leading edge value of x at the wing root to .5 semi spans 
ahead of the leading edge in steps of .05 semispans, and z varying from .25 
semispans above and below the leading edge z coordinate. The upwash is calcu- 
ated using all 18 trailing vortices without numerical difficulty. It should 
be pointed out however that the upwash values at the leading edge are exces- 
sively large because of the proximity to the bound vortex at the wing quarter- 
chord . 

The WASH program is made up of a mainline, four subroutines, and a 
function subroutine. The mainline is the control center for reading and 
writing the specified and calculated information as well as the coding 
responsible for tracing the position of the vortex system downstream. The 
mainline calls the following subroutines: (1) SPLINE to obtain a spline fit 

of the bound vortex strength versus semispan distance, (2) WB which computes 
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the induced velocity components at a point due to the bound vortex using the 
function WBB, (3) WTRAIL to calculate the induced velocity components at a 
point due to a stra i.ght— I i ne trailing vortex segment, and (4) PLANE to 
determine the downwash information in the plane specified by the user. The 
mainline contains the coding responsible for tracing the shed vortices down- 
stream using a maximum of 18 steps behind the trailing edge. Using this 
deflected vortex system the downwash information is calculated at each of the 
grid points given by subroutine PLANE for the user specified plane of interest. 

The program listing, sample data set, and sample output are presented 
along with User Instructions in Appendix A. The next section of this report 
presents a brief discussion of results obtained by running the program for 
several test cases, some of which are compared with experimental data. 
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DISCUSSION OF PROGRAM RESULTS 


Results given by the WASH program were compared with two separate examples 
of USA 45 wings (Ref. 7) at different angles of attack, and, in addition, a 
sample calculation of the wake sheet location and downwash at the' tail was 
obtained for the ATLIT (Advanced Technology Light Twin) aircraft. The theo- 
retical section lift distribution taken from Reference 7 and given in Figure 
16 served as the basis of the calculations in the WASH program. Figures 17 
and 2 show the experimental lines of constant downwash angle in the plane of 
symmetry behind the wing and the angles predicted in Reference 7, respectively. 
The results of the WASH program are given in Figure 18. Reference 7, notes that 
the theoretical downwash of Figure 2 along the wake centerline at large dis- 
tances (two semi spans) appears to be approaching a value of 7.8°, while the 
experimental value is less than 6.6°. Reference 7 indicates that the theoret- 
ical value given in Figure 2 is too large because the method did not include 
spanwise curvature. The WASH program, which includes spanwise curvature, 
gives a value closer to the experimental result in the vicinity of the hori- 
zontal tail location. Near the wing trailing edge the theoretical method of 
Reference 7 over-pred i cts the downwash angles while the WASH program under- 
predicts it as compared with the experimental value. It should be noted that 

some difference between theory and experiment may be attributed to the fact 

that the r distribution is theoretical rather than experimental. 

The experimental T distribution given in Figure 19 for the USA 45 airfoil 

was inserted in the WASH program to compare the spanwise variation in downwash 
with experimental values presented in Figure 20. The WASH program results are 
presented in Figure 21. A comparison of these two figures indicate that the 
total flow angular deviations computed by WASH are smaller in general than the 
experimental values; also. Figure 20 depicts positive sidewash angles both 
above and below the wake centerline while WASH predicts positive below and 
negative above as would be expected by vortex theory. 

The wake sheet is given for each USA 45 r distribution discussed above in 
Figures 22 and 23. Each wake sheet figure was obtained by executing the plot 
data set from the WASH program in the PLOT program of Reference 4. While both 
distributions are for the same airfoil section, an examination of Figures 16 
and 19 gives an indication of how wake sheet shape is affected by the section 
lift d i stri but ion . 

NASA is currently supporting the development of the ATLIT aircraft to 
investigate the potential of new aerodynamic techniques for improving the 
performance of light twin engine aircraft. A three- view of the ATLIT is given 
in Figure 24 and a perspective view is given in Figure 25. A contributer to 
this program, N. C. State University, is currently involved in both a perfor- 
mance prediction and a flight evaluation program for the ATLIT, the results 
of which should be published as a NASA Contractor’s Report early in 1977. In 
conjunction with this program, it was decided to analyze the ATLIT wing using 
the WASH program to find the downwash characteristics at the horizontal tail 
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in the plane of symmetry. The section lift coefficient distribution of the 
wing with fuselage was obtained by the methods given in Reference 4. The 
nacelles were neglected for this test case. The program input data set is 
given in Figure A-2. Figure A-3 gives typical program output: tabulated data 

for the spanwise variation of the downwash and sidewash in the plane of the 
horizontal tail, the downwash in the plane of symmetry, and the upwash aheaa 
of the wing for a wing lift coefficient of 1.22. The wake sheet plot (Figure 
26) indicates that the high wing aspect ratio (ten) yields relatively low 
values of downwash in the vicinity of the horizontal tail. Figure 27 is 
included to show the variation of upwash ahead of the wing in the plane of 
symmetry along with the downwash variation behind the wing. The excessively 
large upwash angles at the leading edge are due to the nearness of the bound 
vortex at the wing quarter-chord. 

The cases presented above, although limited in number, show good agree- 
ment with experimental data. These test cases also indicate the type of 
information yielded by the WASH program; to the extent that interference 
effects from the fuselage and horizontal tail can be neglected, the program 
fulfills the design objectives of defining both the wake sheet location in 
the vicinity of the horizontal tail and the downwash distribution in both 
the spanwise and freestream directions. The program’s capability for com- 
puting upwash ahead of the wing has also been demonstrated. 
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Figure 16. Theoretical T(y) distribution for the 3.66 m USA 
45 tapered wing with = 1.175 (Ref. 7). 



Figure 17. Experimental downwash angle contours in the plane of symmetry 
behind a 3.66 m USA 45 tapered wing with C L = 1.175 (Ref. 7). 
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Figure 18. Downwash angle contours in the plane of symmetry behind a 
3.66m USA 45 tapered wing at C|_ = 1.175 using WASH program 
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Figure 19. Experimental span load distribution for the USA 45 wing 
at C, =1 .35 (Ref. 7). 
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Figure 21. Air flow 7.92m bi 
WASH program witl 




Figure .22. 


Wake sheet end view and orthographic projection of the USA 45 wing at C L 
(T distribution taken from Figure 16). 


= 1.175 
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Figure 23. 


Wake sheet end view and orthographic project! 
(T distribution taken from Figure 19). 


on of the USA 45 wing at C L = 1.35 
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Figure 27. Upwa*sh and downwash angles in the plane of symmetry of the 
ATLIT wing at C = 1.22. 




PROGRAM FOR VELOCITY 


PROFILE CALCULATION 



GENERAL THEORY 


Once the location of the wake sheet centerline has been established using 
the WASH program, a viscous solution behind the wing may be obtained by solving 
the two-dimensional boundary layer equations along wake streamlines defined by 
the shed vortex paths. This same technique was applied to boundary layer 
solutions over arbitrary bodies with good success in Reference 4, and the 
results should be reasonable in this application as long as (1) the spanwise 
flow is small compared to the freestream flow, and (2) the streamline curva- 
ture is not large. These conditions of course are violated as one moves down- 
stream from the wing tip and the vortices begin to roll-up. However, for 
streamlines which pass near the horizontal tail the conditions for the use of 
two-dimensional boundary layer analysis are usually satisfied quite well. 

Thus, this approach should be adequate for the purpose of this study which is 
to determine the magnitude and direction of the flow in the vicinity of the 
horizontal tail. The WASH program permits only nine shed vortices to the 
half-wing; thus, probably no more than the three most inboard vortex traces 
will pass inboard of the tail tip section.* 

As noted in the literature review, the inclusion of curvature in the 
boundary layer equations greatly increases the complexity. In order to 
determine whether or not curvature effects should be included in the analysis 
results from several executions of the WASH program were examined. When 
vortices for typical cases of spanwise lift distributions were traced down- 
stream, it was found that while the inboard vortex streamlines were deflected 
in the vertical and spanwise directions, the curvature of these streamlines 
was small due to the relatively slow rate of change of downwash and sidewash 
angles (see Figures 14 and 26) along a given vortex path. Since the curvature 
of the inboard vortices is small, it should not play a dominant role in a 
boundary layer solution near the horizontal tail. Consequently, the curvature 
terms in the boundary layer equations were neglected, greatly simplifying the 
solution process. 

In order to find the flow velocities along streamlines in the vicinity 
of the horizontal tail including the effects of fluid viscosity, the computer 
program WAKE was written. Given an initial velocity orofile at the airfoil 
trailing edge ++ , the nonsimilar two-dimensional boundary layer equations are 
solved using a finite difference method with even steps ize in the transverse 
direction and unequal stepsize in the streamwise direction. The remainder of 
this section will be concerned with a discussion of the WAKE program theory. 

The boundary layer equations (Equations 11 and 12) are written in a form 
more amenable to numerical solution through use of both the stream function \p 
(defined by u = and v = - ■ which automatically solves Equation 11, and 

a proper scaling of the x and y variables. The scaling variables n and £ 
and the stream function ip were assumed to be of the form 

* Because of typical span loadings, the third vortex out from the fuselage will 
probably be positioned at around 0.4 semi-spans, a location which should be 
further outboard than the tip of the horizontal tail plane on most light 
a i rcraft. 

++The computer program providing this initial velocity profile is given in 
Ref. 4. 
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n = — y and Z = - (33) 

/i0Oc.+x C 

d 

and 

ip(x,y) = U^/l OOc^tx f(£,n) (34) 

= wing chord 
= wing thickness 
= wing section drag coefficient 

= freest ream velocity 

= nondimensiona I distance downstream of the trailing edge 
= nond imens iona I , scaled y coordinate normal to the x-direction 
=. function of £ and n which will appear as the unknown in the 
transformed boundary layer equations. 


where 


Uco 

K 

n 

f U,n) 


Since the continuity equation is identically 
the stream function, it is necessary only to 
by Equation 12 and rewritten, using Equation 


satisfied by the definition of 
solve the momentum equation given 
1 3, as 


u , _3u_ _ ‘ 3^u 

3x V 3y e 3 2 ' 

7 y 


(35) 


Introducing the definition of 
scaled variables. Equation 35 


the stream function and replacing x and y by the 
becomes 
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03 


100c .tx 

d 


f ! ! 1 


(36) 


where the primed quantities denote derivatives with respect to the variable n. 
Simplifying the above equation, one has 


11 111 _ + 11 l )f „ = _ e_ . . 

c 3? 2x 3? c U 00 ( 1 0Oc^tx) 


( 37 ) 


or. 


multiplying by x and 


rearranging. 
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e 


( 38 ) 


100c ,+U 

d c 


+ — ff" - £-{ f t^JLL _ f»— > = 
2 rT 85 T 85 ^ 


Equation 38 must therefore be solved at every 5 station (non-simi iar solution) 
downstream in order to obtain the velocity as a function of n. It is a third 
order partial differential equation in f requiring the three boundary 
conditions: 


u(x,y-*°°) 
u(x,y-*— °°) 
v(x,y=0) 


U -* u( 5 ,°°) = U — f ’(5,°°) = 1 

- u ( 5 , -°° ) = — ► f , (5,-°°) = 1 

0 f(5,0) + 25 If = 0 . 

5=0 


(39) 


Equation 38 together with the boundary conditions given in Equations 39 
is solved in the WAKE program by the finite difference technique described 
below. Let f' = T; then 


T" + — T f - Cl 5(f - T»||} = 0 , (40) 

where c. = 100c_,tU /e. 

1 d 00 

Simplify Equation 40 to read 

T" + o^T' + a 2 H = 0 ' (41) 

where 

a 1 = c^ + 5 -|^-) and = - c 1 5f ' . (42) 

Equation 41 has, through this process, been converted into a second order 
"linear" equation in T with supposedly known coefficients. This equation is 
solved by assuming a solution to evaluate the a T s and then, using an iterative 
process, T is computed numerically and the a's are corrected after each 
iteration until convergence is attained. Consider the grid system shown in 
Figure 28. If the solution is known at some station m for all n, then using 
central and backwards difference formulas , the solution is constructed at the 
station m + 1 for ail n. The derivatives are approximated with the following 
relations taken from Reference 19: 
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Figure 28. Grid system for finite difference solution. 
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2 1 • 2 ' 


where ~ ? m _i> h 2 “ ? m+1 ~ and An = equal steps ize in the n direc- 

tion. Equation 41 is now rewritten as 


* This method is based on that presented by Cebeci and Smith in McDonnell- 
Dougias Aircraft Co. Inc. Report No. DAC-67130, October 1968, "A Finite- 
Difference Solution to the Incompressible Turbulent Boundary Layer Equations 
by an Eddy-Viscosity Concept". 
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1 a 1 

c— 1 — + X:t 


-.+ ^ + a,{ u - } JT _ + C— 1— . 'jt 


(in) 2 2An n+1 ,m+1 " (An) 2 "2‘h 2 (h, + h 2 > JJ n,m+1 (An) 2 2An J, n-1 ,m+1 


h 9 h + h 7 

+ «9tKr-T h" Tr i " r -~- }T D =? 0 

2 h ^ C h 1 + h 2 ) n,m-1 n,m 


(44) 


or. 


where 


A T •_ L1 ! , + B T . , + C T . . = D , 

n n+1,m+1 n n,m+1 n n-1,m+1 n 




(45) 
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n >2 a 2 h 0 ( h , + h„) 
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2 1 


( Ar)) 
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n 2 h i^ 1 + h 2 ) n,m-1 h 1 h 2 n,m 


(46) 


The total number, n , of n-stations must be odd. The bottom of the wake 
(rr>-°°) is defined by n = 1 while the top (r}-H-°°) is defined by n = . The 

wake center is defined by n m - . = (n ma;< -■ 1)/2 + 1. The boundary conditions 
can therefore be specified us i rig the index n as 


n ^ ^n,m+1 ^ 


n = n . , 
mi d 


f(5,0) + 25 


3f 

S5 


= 0 


(47) 


(5,0) 


n = n T . . = 1 

max n,mtl 
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Equation 45, when written in matrix form, becomes; 
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- . - 
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0 
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^2 

A 2 

0 




T 2 


D 2 
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B 3 

A 3 




T 3 


°3 





C 0 B . 

n-2 n-2 

A „ 0 
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T n-2 


°n-2 





0 C . 

n-1 

B i A - 

n-1 n-1 


T n -, 


D , 
n-1 

- 




0 0 

0 1 


T n 


1 

The 

solution 

procedure 

for this set of 

equations can be greatly 

simp 1 if ied 


because the matrix is tri diagonal. If the matrix is denoted by A, then 

A T - D . (49) 


Now let 


where 


A = L X 


(50) 


e 2 
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0 0 
“ 2 0 



3 < co . 0 

n-1 n-1 


0 3_ w 

n n 


0 1 -E 


n-1 


0 0 1 


and 
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0 


0 


0 


LX = 


thus, 


- U) E, 


C — E i ^2 + w 2^ 


e 2 c 2 


* E 2“2 


^ -E 2^3 + w 3^ 


co ^ = 1 


E 3 U 3 


(-E 3 0 4 + u> 4 ) - w 4 E 4 


- to ^ E ^ = 0 


$ 3 - C 3 


- £,62 + « 2 = b 2 


- M2 e 2 = a 2 


b 4 - c 4 


- E 2 $ 3 + oj 3 = B 3 


' “3 E 3 = A 3 


3 = C 

n n 


- E . 0 + to = B 

n-1 n n n 


- to E = A 
n n n 


Applying the above equations, a step by step procedure is used to evaluate the 
0’s, to' s, and E's given the A n ’s, EMs, and C^'s. From Equations 49 and 50 ; 


L X T = D , 


(51) 


and i f 


X T = e , 


then 


L e = D , 


(52) 


(53) 


or 
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If e -j - 
below 


1 and Ej = 0 then all the e ? s and E's are evaluated using the relations 
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(54) 


The matrix Equation 52 may now be written as 
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Starting with T , the n values of T are obtained by using the recursive rela- 
tion given below: 


T n ' 


+ E T 

Vn+1 


( 55 ) 


These T^ values represent the solution at the £ station m + 1 
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The general procedure given above may be summarized in the following 
manner: 

(1) The solution is assumed to be known at the station £ for al I n. 

(2) An approximate solution, T^ ' s assumed at the station in 

order to evaluate the A T s, B ’s, C ' s, D ’ s, e T s, and E ’s. 

n n * n n n n 

(This approximate solution at station m+1 is usually taken to be 

the solution found for station m.) 

(3) Usi.ng the computed values of e and E , the solution T is 

calculated using Equation 55. n n 

(4) T . , and T are then compared to see if they are the same at 

old new r 

every n station to within a certain accuracy; if not the T^^ is 
set equal to T^ ew and steps (2) and (3) are repeated. 

(5) This iteration procedure is continued until T . , and T are 

K old new 

sufficiently close signaling a converged solution at the £ station 
m+1. Usually, no more than two or three iterations are required. 

(6) The procedure given in steps (1) through (5) is merely repeated 
for each subsequent downstream station. 

As noted in (1) above, the finite difference technique requires that a 
solution be known at some £ station in order to compute solutions downstream. 
While an actual velocity profile in the wing wake is not known in general, a 
profile shape can be found at the trailing edge of an airfoil on both the upper 
and lower surfaces whenever the flow over the airfoil is calculated. One 
method of finding the viscous solution over an airfoil is given in the 2-D 
Airfoil program of Reference 4. The flow is computed by solving the boundary 
layer flow over the airfoil with a momentum integral technique. Using the 
program output, the velocity profiles on both surfaces may be determined from 
their respective boundary layer thicknesses and form factors using the equations 
below: 


U 


x_ 

6 


u 


(H u -1)/2 


and 



U 


co 



(H^-D/2 


(56) 


where 6 = boundary layer thickness at the trailing edge 

H = turbulent form factor of the trailing edge 
u & l - subscripts denoting upper and lower surfaces. 

Although Equations 56 do yield an initial velocity profile, there is still 
a problem due to the singular nature of the boundary layer equations at the 
airfoil trailing edge. From Equation 33, it is evident that n is unbounded at 
£ = 0. Some procedure or technique must therefore be used to provide an initial 
profile downstream of the trailing edge. An investigation of the experimental 
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data given in Reference 26 and 27 indicated that at small distances downstream 
of the trailing edge, the velocity profile shape is still much the same as at 
the trailing edge. Based on this information, the assumption was made that the 
profile shape at a small distance downstream, say .01 chord lengths, .is the 
same as the profile at the trailing edge. This technique was used to establish 
a profile at £ = .01 in the WAKE program, and by assuming this profile satis- 
fies the wake boundary layer equations, the initial requirements for the finite 
difference solution are satisfied. Although the program permits initial values 
other than .01 by requiring the value to be read into the program, this value 
was found to be satisfactory for all test cases presented in a later section of 
this report. Values of £ less than .01 tended to make the maximum value of n 
too large, requiring an excessive number of steps in the n direction (Equation 
33). For values larger than .01 the profile shape can no longer be approxi- 
mated adequately by the trailing edge profile shape as seen from the experi- 
menta I prof i I es. 

Another problem which must be discussed arises from the fact that the 
boundary layer equation (Equation 35) contains no terms to account for pressure 
fluctuations in the y-direction. In the development of the boundary layer 
equations the pressure was assumed constant in the direction normal to the flow, 
eliminating the pressure gradient term. Far downstream this pressure variation 
is obviously small. Near the airfoil trailing edge, however, the deviations 
from freestream static pressure may be large. The fluctuations are virtually 
impossible to predict unless the complete Navier Stokes equations are solved 
which at the present time is an unrealistic task. The WAKE program was written 
with the pressure terms excluded with the hope that it would still match exper- 
imental data reasonably well. The results obtained indicated that these hopes 
were attained. Correlations with several test cases are given in the next 
section of this report. 

The eddy viscosity model used in the WAKE program is that given by 

Equation 13 in the literature review. This model was derived by Prandtl for a 

fully-developed wake flow . 'In general, the velocity term in parentheses 

represents the difference between the maximum and minimum velocity in a profile 

while b represents the wake thickness or height at a velocity station halfway 

between the maximum and minimum velocity values. These parameters can be 

estimated given a velocity profile in a fully-developed wake, but if the 

initial profile is the trailing edge profile, then another set of parameters 

must be defined which will correlate with experimental measurement. It should 

be emphasized that while some experimental data exist for velocity profiles in 

the wakes of airfoils, these data must be described as very limited (Refs. 7, 

26, 27, 28). Noting Figure 29, the velocity difference for the trailing edge 

profile would always be u since u . is zero, while in the fully-developed 
* ' mgx min w / 1 

profile it is never as large as u ma ^,. Also, the wake width at half-depth for 

the trailing edge profile is always smaller than the same width for a fully^- 

developed wake. 
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Figure 29. Sample velocity profiles at airfoil trailing edge and in 
a fully developed wake. 


A value for the constant, c, in the eddy viscosity model must also be 
chosen which will bring about a match between experimental and theoretical 
data. In Equation 40 the eddy viscosity appears in the parameter c^ as 


c. = 100c, tU /e where e = bc(u - u . ) . 

1 d oo max min 


This may also be written in the following form 


c. = 100c, t/e where e = bc(u - u . )/U 

1 d max mm 00 


(57) 


For wakes considered in this report, the maximum velocity in the profile will 

u m i n 

be the outer edge velocity, U^. Then e = bc( 1 - — -) . If c is defined as a 

positive constant, then the value of c^ is obviously positive. By making 

several test runs of the WAKE program, it was found that at a given distance 

downstream of the initial profile the larger the value of c^, the lower the 

value of the centerline velocity and the smaller the wake width. Thus, it was 

necessary to find new definitions of the terms b and ( 1 - u . /U ) to be used 

m 1 n 00 
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with the trailing edge profile shape to evaluate e. New definitions were 
sought which would, in general, produce the proper variation of e with wing 
angle of attack without making the functional relationships more complex. As 
seen in Figure 29, in order to achieve approximately the same wake widths 
using these two profiles, the width of the trailing edge profile must be 
measured at a y station with a larger velocity. Using this new wake width, 
a relationship was sought between the term ( 1 -u . /U) and the profile drag 
coefficient. Based on a comparison of many computer runs of the WAKE program 
and their correlations with experimental profiles, it was found that (1) if 
b is defined as the wake width at the y station where the velocity is 75 
percent of the freestream value, and (2) if u . /U is represented by the 
expression (1 - c^ / c ^1, then good agreement^nth^downstream wake profiles 

is obtained. min 

Having defined new relationships for b and (1 - u . /U ), the only remain- 
ing parameter needed to complete the definition of e is the constant c. While 
the value of c for wakes behind circular cylinders using the conventional 
definitions of b and (u - u . ) is qiven as .047 in Reference 20, little 
has been published in tne open literature regarding the value of c in airfoil 
wakes. For profiles generated using the WASH program with the new parameters 
defining the eddy viscosity, a constant value of c = .03 was found to work 
very well for all the airfoils tested. Although this value produces good 
correlation, the value of c is read into the WAKE program to allow the user 
the option of modifying the constant if -he desires. The parameters b and 
(1 - c^ / c are also defined on cards FIN 51 and 52, respectively, in the 

program m I Vst i ng given in Appendix B. This permits the user to experiment 
with new ideas when more experimental data become available. 

it is well to note here that a constant value of e for all x locations 
is realistic only for relatively weak wakes, -> 0. e is really a 

measure of the eddy size and intensity in the wake. In high drag situations, 
large eddies are shed. These have very large internal shears which cause 
them to disintegrate progressively into smaller eddies with smaller shears 
as one moves downstream. For such a situation a constant eddy viscosity is 
obviously not realistic. It will be noted, however, that the procedure for 
computing the boundary layer on the wing does not admit significant regions 
of separated flow. Since the wing drag is always relatively small when the 
boundary layers are attached and since the airfoil boundary layers are com- 
puted with constant eddy viscosities, it seems reasonable to continue this 
practice in the wake computations. As a result, however, the determination 
of the momentum defect must be limited to those cases where these conditions 
of unseparated flow and low drag are always satisfied: generally for section 

lift coefficients of 0.8 or less. 

In order to carry out the numerical integration of the WASH program, the 
user must specify the grid spacing in both the £ and n direction. The program 
is designed to use a maximum of 2001 points in the n direction or 1000 points 
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on each side of the centerline. The maximum values of n for the initial upper 
and lower profiles are computed using the boundary layer thickness 6- in 
Equation 33: 


(n 


) 

max 


upper 

or 

lower 


5 

upper 

or 

I ower 


1 00c .tx 
d 


initial 


(58) 


Dividing the larger of the two n values by 1000 yields the minimum stepsize 
in n. While the program is written for the general case (unequal stepsize) 
it is recommended that only even stepsize be used. If the specified even 
stepsize is too large, very small oscillations in u/ll^ at the outer edge of 
the wake profile may result; hence, a stepsize which increases with n would 
only compound the oscillations. For every case presented in this report a 
stepsize of .00125 was found to be satisfactory. Using this stepsize and n 
from Equation 58, the number of points in the n direction may be computed froft 


NPTS = 2k. n /. 00125 + 1 
1 max 


where NPTS is the total number of n stations for the entire wake and k. is a 

factor to insure that an n region is defined which is a little larger Than 

the actual wake thickness at the initial station; k^ should therefore have a 
va I ue of 1.1 or 1.2. NPTS must be an odd integer. For the cases presented 
in this report it was usually specified as 1201, 1401, 1601, 1801, or 2001, 
depending on the maximum value of n. Since the value of NPTS is directly 

proportional to execution time it is suggested that less than 2001 points be 

used if possible; however, inexperienced users will probably want to choose 
one of the suggested values of NPTS given above depending on the value of 
n 

max 


Along with an n stepsize, a stepsize must also be specified in the 
streamwise or K direction. For 5 an uneven stepsize which is very small 
initially and progressively larger downstream is recommended. The WAKE 
program is designed to use a stepsize which increases every step by the factor 
(1 + DXSTEP) . If DXSTEP = 0 then equal step downstream will be taken. Unequal 
steps are suggested because the initial solutions are very sensitive to step- 
size while downstream they are relatively insensitive. The user must specify 
both the initial stepsize DX and the stepsize increment DXSTEP. For most of 
the cases presented here, an initial stepsize in 5 of .00005 was used with 
DXSTEP = .05; thus, at every step the stepsize increased by 5 percent from its 
previous value. Using this stepsize and increment, approximately 150 steps 
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are required to traverse a distance of two chord lengths downstream (£ = 2.0). 
One disadvantage of the unequal stepsize is that- if the solution is desired 
at £ = 2.0 it is very difficult to choose DX and DXSTEP so that one of the 
steps falls on this £ station; the solution nearest the desired station must 
therefore be used. Figure 30 presents a list of downstream station number 
(NXPRNT in the WAKE program) and £ distance for the different stepsize 
increments .05, ,075, and .10 all starting with an initial stepsize of .00005. 
This data should be helpful in determining the station number at which a 
solution is printed as well as the total number of steps required downstream 
(NXSTEP in the program) to achieve a desired £ distance. 

The WAKE program consists of a mainline and five subroutines. The main- 
line (1) reads and prints the input data set, (2) calls subroutine PR0FIL to 
obtain the initial profile at £ = .01, (3) calls subroutine GETF to obtain a 
value of f (Equation 41) at the first £ station downstream given the derivative 
f T or T, (4) calls subroutines ABCR, C0EFF, and SOLVE to obtain a new value 
for the profile T at the £ station, (5) compares the new profile with the last 
profile calculated at the £ station to see if convergence is attained, and 
(6) after convergence, prints the solution and repeats the same procedure for 
the next downstream station. A listing of the WAKE program with user instruc- 
tions, a sample data set, and sample output are given in Appendix B. The 
sample data set should be particularly valuable to the unfamiliar user in 
giving examples of stepsizes and possible print options. 
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DISTANCE FROM TRAILING EDGE IN CHCRDS 

ST AT ION 

NUMBER DXSTEPx.OS OXSTEP* *0 75 OXSTEP* ■ 0 l 


DISTANCE FROM TRAILING EDGE IN CHORDS 

station 

NUMBER DxSTEP=.05 OXSTEP*-075 DxSTEP*.Ot 


1 

0.010000 

0.010000 

0.010000 

84 

0.066374 

0.278990 

1.37271 0 

2 

0.01 0050 

0.0 10050 

0.010050 

85 

0 • 069242 

0*299214 

1.509031 

3 

0.010103 

0.010104 

0.010 105 

86 

0.072254 

0.320955 

1 .658985 

4 

0.010158 

0 .0 10 162 

0.010166 

87 

0*075417 

0*344327 

1.823933 

5 

0.010216 

0*010224 

0.010232 

08 

0.078738 

0.369452 

2.005376 

6 

0.01 0276 

0.010290 

0.010305 

89 

0.082225 

0*396460 

2.204964 

7 

0.010340 

0.0 10362 

0.010386 

90 

0*085886 

0*425495 

2.424910 

S 

0.010407 

0.010*39 

0.0 10*7 * 

91 

0* 089730 

0.456707 

2.66601 l 

9 

0.010477 

0.0 10522 

0.010572 

92 

0. 093767 

0*490260 

2*931662 

to 

0.010551 

0.01061 t 

0.010679 

93 

0 • 098005 

0*526330 

3.223879 

1 t 

0.010629 

0.010707 

0.010797 

94 

0.102455 

0.565104 

3.545317 

12 

0.010710 

0.010810 

0.0 10927 

95 

0.107128 

0*606707 

3.898898 

13 

0.010796 

0*0 10921 

0.011 069 

96 

0.11203S 

0*651596 

4. 287038 

14 

0.010806 

0.011040 

0.011226 

97 

0.11 7186 

0.699766 

4*715672 

15 

0.010900 

0.011168 

0.011399 

98 

0.122596 

0.751548 

5. 186289 

16 

0.011079 

0.011306 

0.011589 

99 

0* 128276 

0 • 807 2 1 S 

5.70346a 

l 7 

0.01 1 103 

0.0 11454 

0.01 1 797 

100 

0. 134239 

0.867056 

6.273415 

1 a 

0.01 1292 

0.011613 

0.012027 

10 1 

0.140501 

0.931305 

6*899006 

1 9 

0.011*07 

0.01 1 70* 

0.012280 

102 

0*147076 

1 .000539 

7.586837 

20 

0.01 1527 

0.011968 

0.012558 

103 

0.153980 

1 .074879 

0.346770 

21 

0.01 1653 

0 .0 12 165 

0.0 12064 

104 

0. 161229 

1.1 54795 

9. 1 60498 

22 

0.011706 

0.012378 

0.013200 

105 

0. 168841 

1.240705 

10.097597 

23 

0 .011925 

0 .Q 1 2606 

0.013570 

106 

0. 1 76833 

1.333057 

1 1.106407 

24 

0.012072 

0.012051 

0.013977 

1 07 

0 • 185224 

1 .432337 

12.216098 

25 

0.012225 

0.0131 15 

0.01**25 

1 06 

0.194035 

1 .539062 

1 3.436757 

26 

0.012306 

0.013399 

0.01*917 

109 

0*203287 

1.653792 

14*779403 

27 

0.01 2556 

0.01370* 

0.015459 

110 

0*213002 

1.777126 

16.256482 

20 

0.012733 

0.0 140J2 

0.016055 

11 1 

0. 223202 

1.909711 

1 7. 081 180 

29 

0.01 2970 

0 .0 14304 

0.016710 

112 

0.233912 

2.052239 

19.668348 

30 

0.0131 16 

0.0 14 763 

0.017432 

113 

0.245157 

2.205457 

21.634232 

31 

0.013322 

0.015170 

0.018225 

1 1 4 

0.256965 

2.370166 

23.796706 

32 

0.013538 

0.015008 

0.019097 

115 

0. 269363 

2.547228 

26. 1 75426 

33 

0.013765 

0.016078 

0.020057 

1 1 6 

0.202382 

2.737571 

28.79201 9 

34 

0.01*003 

0.0 16584 

0.021 1 1 3 

117 

0. 296051 

2.9421 88 

31 .670271 

35 

0.014253 

0.017128 

0.022274 

1 18 

0.310403 

3. 1 62152 

34.036348 

36 

0.01451 6 

0.0 1 7713 

0.023551 

(19 

0.325473 

3.390614 

38. 319033 

37 

0.014792 

0.0 10341 

0.024956 

120 

0.341297 

3 .652010 

42. 149906 

30 

0.015001 

0 .0 190 1 7 

0.026502 

12 l 

0.357912 

3.926071 

46.364034 

39 

0.015305 

0.0 1974 3 

0.028202 

122 

0.375358 

4 . 2 t 9026 

50.999400 

40 

0.015705 

0.020524 

0.030072 

123 

0.393675 

4.535613 

56.090407 

41 

0.0100*0 

0.021303 

0.032 1 30 

124 

0.412909 

4.875084 

61. 707385 

42 

0.016392 

0.022265 

0.034393 

125 

0.43310S 

5.24001 5 

67.077174 

43 

0.010767 

0.023235 

0.036882 

126 

0.454310 

5.63231 6 

74.663941 

44 

0.017150 

0.024278 

0.039620 

127 

0.476575 

6.054040 

82* 1 29385 

♦ 5 

0-01 7557 

0.025 396 

0.042632 

128 

0 . 499954 

6.507393 

90.341374 

♦ 6 

0.01 7905 

0.026603 

0.045945 

129 

0. 524502 

6.994748 

99. 374561 

47 

0.018434 

0.027898 

0.049590 

130 

0.550277 

7.510654 


48 

0.018906 

0.029291 

0.053599 

131 

0.577341 

0*061853 


49 

0.019401 

0 .0 30788 

0.058009 

132 

0.605758 

0.607292 


50 

0. 019921 

0.032 397 

0.062959 

133 

0.635596 

9. 338 136 


5 1 

0.020467 

0.0 3* 1 26 

0. 068195 

1 34 

0.666926 

10*037799 


52 

0.021041 

0.0 35986 

0.074065 

135 

0.699822 

10*789934 


53 

0.021 643 

0.0 37985 

0.000521 

1 36 

0. 73436J 

1 l .690479 


54 

0.022275 

0.040 1 34 

0 .087624 

1 37 

0.770631 

12.467665 


55 

0.022939 

0.0 42444 

0.095436 

138 

0.80871 3 

13* 402040 


56 

0 .023636 

0.0*4927 

0.10*030 

1 39 

0.848698 

14.406493 


57 

0 .07*307 

0.047597 

0. 1 13403 

140 

0.890683 

1 5.406279 


50 

0.0251 36 

0.0 50466 

0 • 123801 

14 1 

0.934767 

16.647050 


59 

0.025943 

0.0 5355 1 

0. 135319 

1 42 

0.901056 

17.094879 


60 

0.026790 

0.056868 

0. 1 4790 1 

143 

1 * 029659 

19.236295 


61 

0.027679 

0.060433 

0.161741 

144 

1 .000691 

20.67031 7 


62 

0.028613 

0.064265 

0. 176965 

145 

1.134276 

22*228491 


63 

0.029594 

0.068385 

0. 1 9371 1 

146 

1 . 190540 

23.094920 


64 

0 .030623 

0.072614 

0.212133 

147 

1 . 24961 7 

25.606348 


65 

0.031705 

0.077575 

0.232396 

148 

1 .31 1640 

27*6121 24 


66 

0.032840 

0 • 0 8269 3 

0.254685 

149 

1 . 376700 

29*662333 


67 

0.034032 

0.000195 

0.279204 

150 

1 .445169 

31 *907008 


68 

0.035203 

0.094 | 10 

0. 30617* 

151 

1.51 6977 

34.300193 


69 

0.036590 

0. k 00468 

0 • 335842 

152 

1 . 59237 . 

36.072008 


70 

0.037978 

0.107303 

0.368476 

1 S3 

1 .671545 

39.6 36709 


71 

0.039426 

0. 1 14651 

0.404373 

154 

1 . 754672 

42.600762 


72 

0.040940 

0 . 1 22550 

0.443061 

155 

1 .841956 

*5.803719 


73 

0 .0*25*5 

0. 13104 1 

0. * 5729 7 

1 56 

1 .93360* 

49.238290 


74 

0.044222 

0. 140 169 

0.535077 

157 

2.029034 

52*930470 


75 

0.045984 

0 .149962 

0 .507634 

158 

2. 1 30876 

56.099555 


76 

0.04 7033 

0. 1 60530 

0*645440 

159 

2*236970 

61 * 166322 


77 

0.049774 

0.1 71870 

0.7090*2 

160 

2* 340368 

65.753096 


78 

0.051B13 

0.184061 

0.778997 

161 

2.465336 

70.663676 


79 

0.053954 

0. 1 97165 

0.855946 

162 

2.5881 S3 

75.984469 


80 

0.056201 

0.211 252 

0.940591 

163 

2. 71 71 1 1 

8K .602604 


8t 

0.058561 

0.2 26 396 

1 .033700 

164 

2.852516 

87. 808100 


82 

0.061040 

0.242676 

1.136120 

165 

2.994692 

94.393007 


83 

0.063641 

0.260 177 

1.248782 

166 

3*143977 




Figure 30. 
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DISCUSSION OF PROGRAM RESULTS 


Correlations of the WAKE program results with available experimental 
velocity profiles are presented in Figures 31 through 37. The figures show 
the experimental data, results obtained by the theoretical method of Reference 
7, and the WAKE program results for specified stations downstream of the 
trailing edge. Each airfoil was analyzed using an initial x/c value of .01 
and an eddy viscosity constant of .03. The profiles presented are not intended 
to show the relative distance of the profile above or below the trailing edge, 
but only the profile shape about the wake centerline as it moves downstream.. 

In Figures 31 through 34 the WAKE program results are shown for two Reynolds 
numbers. In order to achieve a theoretical drag coefficient approximately 
equal to the experimental drag measured in References 26 and 27, it was neces- 
sary to run the two-dimensional airfoil program of Reference 4 at a Reynolds 
number of approximate ly half (.21 million) the indicated test value of .42 
million. Tunnel turbulence or model surface conditions during the wind tunnel 
test may have produced more drag on the airfoil than would normally be expected 
at the test Reynolds number. The experimental drag was measured using a 
momentum rake which indicates the momentum defect in the wake at a station 
downstream of the trailing edge. Since the velocity profiles depend on drag 
coefficient, the calculated profiles should agree better with experiment if a 
theoretical drag coefficient is used which agrees with the experimental value. 
It was therefore decided to compute and compare profiles using both the 
specified test Reynolds number and a "corrected" Reynolds number, _i_.e^ , one 
at which the theoretical drag matches the experimental drag. 

Figures 31 and 32 show that the WAKE program results agree well with 
experiment for the symmetrical Joukowski airfoil at both angles of attack (0° 
and 6°). While results at both Reynolds numbers show good agreement, that 
for the lower Reynolds number case is better. Although no experimental data 
were available at x/c = 2.0, the results are shown at this station for this 
and all other airfoils to give the reader an idea of the profile shape in the 
vicinity of the horizontal tail. 

Figures 33 and 34 present program results for a symmetrical Piercy air- 
foil at angles of attack of 0° and 6°, respectively. The experimental points 
were taken from Reference 27. The WAKE program profiles give good agreement 
with experiment except near the minimum velocity region of the wake centerline. 
While the method of Reference 7 predicts the minimum velocity quite well, it 
over-pred i cts the wake thickness at other velocity stations in the profiles. 

Experimental measurements of a modified 15 percent symmetrical airfoil 
at an angle of attack of 8° are presented in Reference 28, and these results 
are compared with the WAKE program prediction in Figure 35. Both the WAKE 
program and the method of Reference 7 predict the profile shape accurately as 
compared with experimental data except near the outer edge where the experi- 
mental data points appear questionable. 
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As a final test case the USA 45 airfoil was analyzed at angles of attack 
of - 5.3° and 1.6° and compared with experimental data presented in Reference 
7. It should be pointed out that the profiles given in Reference 7 were not 
presented to make comparisons of this type, but rather only to give approxi- 
mate dynamic pressure profile shapes in the general sense. The reader should 
bear this in mind when considering these two figures because the experimental 
profile shapes in Reference 7 have been greatly enlarged to make these com- 
parisons. For the angle of attack - 5.3° the WAKE program did a better job 
of predicting the profile shape; however, for a = 1.6° the method of Reference 
7 appears to match experiment better. 

In general, agreement between the WAKE program profiles and experiment 
was good and, in most cases, better than the method of Reference 7. The 
minimum velocity at x/c stations near the trailing edge was generally larger 
than experiment while centerline velocities at stations greater than 1.0 were 
generally smaller than experimental values. The correlation of profile width 
using the WAKE program with experimental profile shape was excellent. The 
reasons the coefficients of the eddy viscosity were chosen to give a better 
fit to the experimental data further downstream of the wing rather than at the 
trailing edge were first that downstream is the region of major interest and 
second that a boundary- I ayer- type analysis obviously cannot describe the 
initial portion of a base flow wake adequately. Thus to favor a data match 
in a region where the analytical model is known to be inadequate is to lose 
the opportunity to attain a match in a region where it should be reasonably 
correct. 
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Figure 35. Comparison of velocity profiles with experiment (Ref. 28) for a modified 
0015 airfoi I at a = 8°. 
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APPLICATION TO THE CALCULATION OF HORIZONTAL TAIL ONSET FLOW 


in the preceding, the programs for determining the wake sheet location 
and the velocity profile in the wake have been discussed independently; con- 
sideration will now be given to combining them so as to compute the onset flow 
field at an aircraft’s horizontal tail. The trailing vortex streamline 
locations are calculated using WASH. Naturally velocity profile solutions 
need be computed only along those streamlines which pass near the horizontal 
tail. Given the section drag coefficient of the wing at the spanwise origin 
of each streamline used in the computation and the total streamline length to 
the x-station of interest, WAKE determines the velocity profile centered about 
the z-location of the wake sheet. The total distance traced by a streamline 
in going from the y-axis to some downstream x-station is given in the WASH 
program output in the section specifying the coordinates of the trailing 
vortex sheet. An example of how these are properly combined is presented in 
Figure 38 for the ATLIT wing at three stations downstream of the wing trailing 
edge in the wing plane of symmetry. The profile at the horizontal tail is 
therefore computed as the profile at the station x/c = 2.0 of Figure 38. 
Similar profiles may be calculated at other spanwise stations in order to 
show the spanwise variation of the profiles. 

In the strictest sense profiles should be calculated along the 4 or 5 
most inboard vortices to describe the spanwise variation of the onset flow at 
the horizontal tail. However, from a practical point of view unless the span- 
wise curvature of the wake is large there is, as may be seen in the wake sheet 
plots, little vertical variation in sheet location for v distance equal to 
normal ta i I semispan lengths . Also, a check of Figures 31 through 37 will 
indicate that at x/c = 2.0 the velocity profiles across the wake are very 
nearly the same for all airfoils tested, indicating that the effect of span- 
wise variations in section drag coefficient on the velocity profiles are 
small. As a result velocity profiles at different spanwise stations will 
usually have essentially the same velocity deficit. If the spanwise and 
downstream curvature is small for vortices which pass in the vicinity of the 
horizontal tail, it is therefore possible to execute only one case of each 
program and obtain the approximate wake position, flow angularity, and profile 
shape in the region of the horizontal tail (the most inboard vortex is used). 
This approximation can be improved if desired by applying the single WAKE 
program velocity profile to each of the desired spanwise vortex positions 
calculated in the WASH program. If a more rigorous analysis is desired, a 
boundary layer solution may be determined along each of the 4 or 5 most 
inboard streamlines noting the spanwise variation of section drag coefficient. 
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Chord lengths above and below wake center line 
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With the profile shape in the tail region determined, estimation of the 
forces on the tail plane due to the flow angularity and velocity variations 
can be undertaken. By executing the WASH program for a range of aircraft lift 
coefficients the wake position as a function of C. may be used to determine 
whether or not the tail ever passes through the wing wake.* If the tail always 
remains 10 to 20 percent of a wing chord, length above or below the wake, 
the velocity deficit analysis may be neglected. Even If the tail is always 
out of the wake, the downwash computations are still important because the 
change in downwash at the tail with angle of attack must be known in order to 
predict tail trim angles and longitudinal stability. If the tail does pass 
through the wake, then the tail must be sized so that it produces the required 
control force even when the toil is located directly in the wake. The effect 
of spanwise variation in downwash may be accounted for by assuming it to be 
similar to a wing twist. While most light aircraft horizontal tails do not 
have twist, a twist variation can be assumed which will produce a spanwise 
angle of attack variation similar to that imposed by the downwash variation. 
This technique can be applied to the three-dimensional program given in 
Reference 4 when computing tail characteristics. 


The determination of the proper flow velocity to use in calculating the 
aerodynamic characteristics of a section of the horizontal tail flying in a 
wake is complicated by the fact that the wake momentum defect is confined to 
a rather thin region which is on the order of the thickness of the airfoil 
for airfoils at low 7 to-moderate angles of attack. In this thin region the 
magnitude of the flow velocity may vary from 0.88 freestream to freestream. 
Since two-dimensional airfoil analysis has not yet progressed to the point of 
calculating aerodynamic coefficients in shear flow, it is necessary to find 
an average velocity squared over the thickness of the airfoil beginning at the 
point on the velocity profile intersected by a line parallel to the airfoil 
chord at t/2. This is best done graphically using plots of the ±pV^ distri- 
bution. 


n 

If the average q = ±pV varies significantly along the tailplane, its 
effects can be included in the three-dimensional aerodynamics calculation 
program (Ref. 4) in several ways: (1) artificially reduce the wing chord in 

proportion to 3 | oca | /Q *, * (2) reduce only the local Reynolds number by the factor 

V. . / V or some function thereof, and (3) multiply C. and C_ by q, . /q . 

local 00 r 1 L D 1 ^ local 

If these procedures result in significant changes in planform or aerodynamic 

properties, they cannot be accommodated readily within the framework of the 

3-D characteristics program. However, if the changes are small enough to be 

capable of linear representation, then one of these procedures will serve 

nicely. Which procedure will yield the best result cannot be stated with 

certainty since time did not permit them to be investigated In detail. The 

author is inclined to favor the first as representing the least compromise 

with physical reality within the context of the two-dimensional airfoil 

analysis. Since the wake region is quite turbulent, it seems reasonable to 

insist that when the tailplane enters it, the analysis should consider the 

airfoil boundary layer fully turbulent from nose to trailing edge. 


* Although the analysis does not include the effects of the tail upwash on the 
wing downwash field, this distortion is not large. The wake position or 
magnitude is therefore not altered significantly. 
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CONCLUDING REMARKS 


The present work provides two computer programs useful for predicting 
the flow magnitude and direction in the region where one would wish to locate 
an aircraft's horizontal tail (pitch trimming surface). The calculation of 
the contribution of this surface te the overall aircraft lift, drag, and 
moment is much enhanced if one has available more realistic values of the 
onset flow. 

WASH ties together for the first time in a computer program a soundly- 
based inviscid method for determining the shape and location of the downwash 
wake sheet with a method previously used to determine the spanwise lift dis- 
tribution for unswept moderate-to-h igh aspect ratio wings and wing bodies in 
incompressible flow. Excellent qualitative agreement is obtained with avail- 
able experimental results. Quantitatively, the predicted flow directions at 
any point in the downwash field agree to within a degree or so of the measured 
values. Considering the experimental error present and the effects of boundary 
layers on the aft fuselage, this is about as good as one could expect. The 
results show there was little spanwise variation in wake position for normal 
horizontal tail semispans. Execution time for the program is 15-20 seconds 
on an IBM 370/165; it is thus a low cost means of mapping the downwash field 
during design analysis. 

WAKE calculates the momentum defect downstream along the wake centerline. 
It uses as input section drag data developed analytically by the aforementioned 
spanwise lift distribution program. It provides velocities that are within 
2 percent of the measured values at all stations for which data were available 
on two wings with different airfoil sections. Computation time for the flow 
along one streamline is 2 min 42 sec. The results indicate that the momentum 
defect due to wing drag is confined to a region no more than 0.43 m high for 
typical light aircraft in the cruise configuration. The region of noticeable 
defect (> 2 percent less than freestream) is only about 0.21 m high. There 
is little spanwise variation in wake velocity profile for typical light air- 
craft horizontal tail semi spans. Thus, unless there is some compelling reason 
to locate the trim surface near the wake centerline, use of WAKE is not 
requ ? red. 

The flight condition for which this type of information is of most 
interest is approach. In this configuration the wake is large and the 
momentum defect significant. While the inviscid downwash calculation as used 
here applies equally well, the viscous treatment suffers from the fact that 
the method used to determine the wing drag cannot accommodate separated flows. 
Further, flows with large separations can not be considered to possess constant 
eddy viscosities as was assumed here. It is apparent, therefore, that if 
horizontal onset flow is to be computed during approach the problems of the 
analytical determination of the correct aerodynamic characteristics of airfoils 
near stall and the characterization of turbulent flow eddy viscosity magnitude 
and decay rate in terms of the size of the initial separation must first be 
so I ved . 
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APPENDICES 



APPENDIX A - WASH PROGRAM 


User Instructions 


The program is written in FORTRAN IV and is designed to run in double 
precision on an IBM 370-165 computer. This program evaluates (1) the location 
of the wake centerline by tracing the path of 18 vortices shed at the wing 
trailing edge, and (2) the total flow angle, sidewash angle, and downwash 
angle behind an unswept wing in any plane perpendicular to the wake centerline 
at a specified distance aft of the wing trailing edge or in the plane of 
symmetry of the wing. Downwash information in a single plane behind the wing 
using 15 downstream steps requires an average execution time of 20 to 25 
seconds. The program requires the specification of the following input data: 


Card Variab le 

Number Name Var i ab I e Descri ption 


1 TITLE 


2 NSTEP 


2 J PUNCH 


3 ALPHA 

3 B2 

3 TR 

3 CR 

3 DX 


4-5 Y 


6-7 CLS 


The 80 character array which is used as a header to identify 
output. Termination of execution is achieved by following 
the last set of wing data to be analyzed by a title card 
having only the word END in the fi rst three spaces. 

The number of downstream steps to be taken aft of the wing 
trailing edge. The total distance the vortex system is 
traced downstream is given by NSTEP*DX. The maximum value 
of NSTEP is 18. 

The control variable giving the user the option of obtaining 
a punched plot data set which may be used to plot the vortex 
system using the PLOT program in Reference 4. JPUNCH=1 
gives punched output while the default JPUNCH=0 gives none. 

The wing angle of attack in degrees. 

The wing semi span, in feet. 

The wing taper ratio (tip chord/ root chord). 

The wing root chord in feet. 

The step size in the streamwise direction specified as ratio 
of wing semispan. 

The 10 variable array specifying the Y-coordi nates of the 
positive Y-axis (axis in the spanwise direction) at which 
the section lift coefficients are read. The order of spec- 
ification is from positive wing tip to wing root as a ratio 
of wing semi span (from Y ( 1 ) at the tip to Y(10) at the root). 
All of the Y values must be non-zero . 

The 10 variable array specifying the section lift coeffi- 
cients correspond i ng to the Y-coord i nates described above 
using the same ordering scheme. 
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Card Vari ab le 


Number Name 


Vari ab I e Descri ption 


8 NXORY 


8 STAT 


Last Card 


The control variable used to specify the type of plane in 
which downwash calculations are made. 

NXORY=0 gives no downwash analysis in a plane. 

NX0RY=1 gives downwash information in a Y-Z plane where X 
is the streamwise direction. The Y-Z plane is located at a 
value of X = STAT (see next variable). Downwash information 
is calculated for the 6 most inboard Y-stations at several 
Z locations. 

NX0RY=2 gives downwash information in the X-Z plane at Y=0 
for X values ranging from .25 to 1.75 semi spans aft. of the 
trailing edge. Z ranges from .5 semispans above and below 
a specified midpoint value of Z = STAT (see next variable). 
NX0RY=3 gives downwash (or negative upwash) information in 
X~Z plane at Y=0 for negative values of X ranging from the. 
wing leading edge to -.5 semi spans ahead of the wing leading 
edge. Z ranges from .25 semispans above to .25 semi spans 
below the leading edge value of Z. 

Variable used with NXORY to specify in semispans the proper 
location of the plane of interest. If NX0RY=1 then STAT 
specifies the X-station at which the Y-Z plane is to be 
located. This distance in semi spans may be either positive 
or negative depending on whether downwash or upwash infor- 
mation is desired, but it must be measured from the wing 
trailing edge. If NX0RY=2 then STAT specifies in semispans 
the midpoint of the Z-range for the downwash calculations in 
the X-Z plane. This Z-range is specified above (positive) 
or below (negative) the trailing edge in semispans. IF 
NX0RY=3 the variable STAT is not used. 

Blank. 


Cards 1 through 8 represent a complete data set for the WASH Program. 
Several data sets may be executed for a s i ng I e program- comp i I e by placing 
subsequent data sets behind the first. More than one plane may be analyzed 
during a program execution by specifying new values for NXORY and STAT on 
other cards placed behind card eight and before the blank card. 

The entire program execution is terminated when the END title card is 
encountered. The format specification for the above data is given in Figure 
A-1. A sample data set is given in Figure A-2, and the output of this data 
set is presented in Figure A-3. 
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Last Card 


» n zs » n n h ji n u m » x 1/ > n u ii ■; u i: is u n u <! so si s 


(BLANK CARD) 

m ii u i) ii is is )) i» » so ?> 17 ;s zi zs. is n a n so i; si 


9 10 II 17 13 II IS IS I' II 19 Z0 Zl 11 Z] ZI ZS ZS 11 ZS Z9 H >1 SI 11 » JS X II X 39 <0 II 1Z 1! II IS IS II U 19 SO SI SZ SI SI SS SS SI SI S9 


CLS(IO) 
1 .39108 


/cLS(l) CLS(2) ■ . . 

.48487 .8103 1.0086 1.1299 1.2084 1.26147 1.3004 1-33085 

I ; S I s . 1 I 9 10 II IZ U II IS IS II II II ZS Zl ZZ ZZ H ZS ZS Zl Zl Z9 M 31 11 11 II IS X SI » 99 lOH IZ I! II IS IS •> II 19 SO SI SZ 51 SI SS X SI 91 S9 SO 61 SZ SJ U SS SS SI 61 69 10 II » « I< 'S IS " I 

/ v/izsl 


' Y(l) Y(2) . . . 

.987946 .9520884 .8933827 .8133888 .7142811 .598922 .4709653 .3355246 

i 7 1 < 5 s 7 I I IQ 11 17 13 H IS 16 If l« 19 ?0 31 27 33 ?i 75 76 2? 2* 29 X Jl J? 33 Jl » 36 JJ 5 39 <0 41 l? U 41 4S 46 4? i| *9 50 SI S7 SJ 5« 55 56 5/ SS ‘•S 10 ll 67 6J M 6S 66 6 ? 61 69 20 71 I? 71 7i TS % ■> It »1 3 


CR DX 

5. 187 .2 

2> 74 JS n 77 ffl 39 70 3: if ■ >: 75 36 37 31 39 *0 4| 47 4J U 4S ii 4' H 49 SO SI S7 Ss M SS 56 SI SS S9 bO 61 67 63 U 55 66 6 ' 6« 69 '0 71 72 71 71 '5 * 'I 71 79 m 


NSTEP J PUNCH 


/^LPHA 

82 T 

/ 10. 1 1 , , 

19.85 

6 9 10 II 1? II 14 IS 16 1’ II n 70 71 


WHITCOMB AIRFOIL CL=1 .2217 B/2 = 19.85 DATA FOR ATLIT 

l ? 3 4 V * 7 « * H) II »? 13 1* IS 16 17 II 1* 20 2! 23 73 71 25 26 77 » 29 M 31 32 33 34 35 36 37 31 3» *0 *M2 II 41 <5 46 47 4 *9 50 Si S? S3 


First Card 


I 4 g 0 1 1 I 0 0 1 1 a I I t M I I I I i I I B It 0 I 1 >M M I I 

1 I I 1 f ■ I 1 1 II 1 1 ■ I 1 I I I 1 1 1 1 1 I I 1 1 1 I 1 M 1 I 1 

722277772777272222222222722222722122: 
3 3 3 1 J3 3 J 1 3 JJ 3 3 3 J 1 3 3 3 3 3 3 3 3 3 3 J 1 3 3 3 3 3 3 3 : 
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5 5 5 5 S 5 5 5 5 5 S 5 5 S 5 5 S 5 S S 5 5 5 5 5 5 S 5 S 5 5 5 S 5 5 5 I 

i I I i M S 6 I 5 i 6 i i i i t S t 8 5 6 6 6 6 6 6 S I 6 6 S t 6 6 6 I 

miimmiiinnmniniminin 

a a 1 1 1 1 a s s 1 1 a a 1 1 a 1 1 1 1 a 1 1 a 1 1 1 1 1 1 a a i a 1 1 
a a a a a a 9 a a a a a a a a a a 9 a a 9 a 9 a 9 9 a a 9 a a s 9 a j 9 


X X ‘0 II IZ I! II IS IS II « IS SO 51 SI 51 51 SS SS SI W 51 60 SI 6? &l 61 16 SS 6! 61 SS lO II 17 11 '< 'S - ■' n SO 

a a a I a o e o o at o o o o o 1 1 1 o i d a 1 1 1 o o i o a o o 9 o o o <i o i o « i 

iMiinninn * 1 1 1 m 1 1 m m 1 n * 1 

2 2 2 2 2 2 2 2 2 2 2 2 7 3 7 7 7 7 7 7 7 2 2 2 2 

3 3 3 3 3 3 3 3 3 3 P i 7 13 3 3 3 3 3 3 3 3 1 

4 4 4 4 4 4 4 4 4 4 /llflHosLiJ fl 4 4 4 4 4 4 4 4 4 4 4 
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55555555555 


SCEEGCS6EE Me HHW 7 66BBEBE6S6BE 

COMPUTING CENTER 

nimnmi \ >. minnniii 

a i a i a 1 1 1 1 1 1 j 1 1 1 1 1 1 1 a a a a i a a a i a i a a a 
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Figure A-2. Example data set for the WASH program 
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Sample Output 


WHITCOMB IIKTOIL CL-I.22I’ 


f)/2«|«,ll DATA FOB ATLIT 


ALPHA 

Ml 


cm 

o* 

tm 

NSHM 

J»ONCH 


O.tOOOOOD 07 
0,1*^000 02 
0.9000000 •• 

o.oioroco (m 

0.7000000 00 
•0.1022010 00 
0.33*0000-01 
IB 


oc cnees 

FEET 

FEET 

SEMI SMtNS 
SEMI SPANS 
SEMI SPAN* 


0.987*460 OO 

0.9920000 oo 

0.09330)0 00 
0.9133090 00 
O.rialMID 00 
0.998*220 00 
0.4709690 00 
0 .33997SO 00 
0.7019000 00 
O.93|9«0O*OI 



O.t 391000 01 


TIB VALUE (V-POSITIVEI 


FOOT VALUE 


ACOAOERCO FULL-SPAN VALUES 


vi r » 

•0.9079400 00 
•0*9920000 00 
-0.0*13830 00 
>0.01 33*90 00 
-0. 71 47010 00 
-0.9969270 00 

-0.3333230 00 
-0.7013800 00 
-0.931 9430-01 



0.9969Z7D 00 
0. 71*2810 00 
0 « 81 33890 00 
0.8933990 00 
0.9920880 00 
0.9879400 00 


0.810)000 00 
O. I 008600 Ot 
0.1129900 Ot 
0.1208400 01 

0 . 13)0 890 Ol 
0. 1 366600 01 
0.1 391080 01 
O. I 391080 01 
0. I 366600 01 
9.1 3 30 8 3D 01 
0. 1300400 01 
0.1261470 01 
0.1208400 01 
0.11 29900 Ot 
0.1006600 01 

0. 4848700 00 


Cl 11 

0.2614640 0 
0.2707280 O 
0.2898990 O 
0.3089610 O 
0 . 3321690 0 
0. 9*1 9690 O 
0.3430260 O 
0.4900170 O 
6.4646220 O 
0.4928330 0 
0.4926330 O 
0.4646220 0 
0.4300170 0 
0.3930260 0 
0.3619690 O 
0. 3321690 O 
0. 9069610 0 
0.2996990 0 
6.2707280 0 
0.2614640 O 


00 


tip vaiue <r - negative! 




if —=>o8i r wet 


LOCATIONS t 


STOEnCTiMS op THE 9 SHEO VOPtlCES ON THE NFGATtVf T-AXIS 
:* muon* I WAGE IS AS9UNFO BOW POSITIVE av IS 1 


-0. 76363SD 
-0.6966020 



GAN5I f 1 
-0.4629740 00 
-0. 9449120 00 
-0.2901900 00 
-0.2790270 00 
-0. 7761 160 00 
-0.7699980 00 
-O. 2929620 00 
-O. 31 3)190 00 
-9. 791 70 10 OO 


TIP VALUE 




BOOT VALUE 


S»L INE PIT Of 



*0.99'' 

-0.900 


-0.660 


UNO V09TPH stre 


0.4 1 7047*«* 
I . 1 I 620746 
1.4191 1899 
1. 60664242 
1 . 77984479 
I .91936986 
2.04247096 
7. 16964] 9 1 
2.76C97997 




0.698 

0.700 

0.790 

0.600 

0.890 

0.900 

1 *080 


7.61 390071 
2. 72137190 
2.67967660 
2.94063701 
3. 0962611 7 
3.17879890 
1.3091 9796 
3.41427966 
3.46349737 
9.90697160 
3.46349737 
3. 41 427999 
9. 1091 9799 
3. I 7679690 
1.096261 1 T 


7.9*471149 
2.700979*7 
2. 16.964133 
7.0*7*7096 
1 .91936989 
1 . 7710*479 
1.60684747 
1.4191 1 099 
I. I 1620246 
0.41 70*749 


Figure A-3. Sample output of WASH program. 
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POSITIVE V-AXIS f F BQX tflNG POOT TO «INC f I n » 


Figure A-3. Continued. 
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LIST IMG OF PLOT DAT* SET FOP SHED VOftTEX STSTCM 


TCOHB AIR 

FH1L CL»I« 

2217 8/2- 

19.65 DAT 

FO* ATLI 

T 






0.2000 

0.4000 

0.6000 

0.6000 

1 .0000 

1 .2000 

1 .*000 

1.8060 

1.6008 

P.8000 

f . 2000 

2. *000 

2. 6000 

2.8000 

3.0000 

3.2000 

3.4000 

3.6900 



A r |c Tt 

0.2696 

0.4032 

0 .53*9 

0 .6566 

0.7638 

0.893* 

0.9227 

0. 9700 



0. 7638 

0. 6566 

0.53*9 

0.4032 

0.2666 

0.147* 

0.0 


e. f» 


0.0 

0.0 

0.0 

6.0 

0 .0 

0. 0 

0.0 


0*0 

D.o 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 


n.-r 

ri , 

0.269* 

0. *0* 2 

0.5360 

0.6300 

0.7659 

0. 0969 

0.9293 

3.97*1 


0-0869 

0.7659 

0.6560 

0.5360 

0.4042 

0.269* 

0- 1*79 

0.0 


-0.0292 

-0.0292 

-0.0233 

-0.0197 

-0.01 73 

-0.013* 

-0.01*0 

-0. 01 26 

-0.6085 


-0. OOPS 

-0.0126 

-0.01*0 

-0.015* 

-0.0173 

-0.0197 

-0.0233 

-8.0292 




0. I*S6 

0.2699 

0. *0*7 

0.5366 

0.6387 

0.7671 

0.69*6 

0.9369 

3.9376 



0.7671 

0 .6587 

0.5366 

0.4047 

0.2699 

0.1*36 




-0.0524 

-0. 0*09 

-0.03*1 

-0.0299 

-0.0267 

-0.02** 

-0.0222 

-9.01*1 


-0.01*1 

-0.0222 

-0. 02** 

-0.0267 

-0.0299 

-0.03*1 

-0.0*09 

-0.052* 

-0.032* 



0. 1 *95 

0 . 270* 

0.4051 

0.5371 

0.6593 

0. ?686 

8.863* 

0.9*87 

3.9337 


n.Mi* 

0.7686 

0.6595 

0.5 371 

.*031 

0.270* 

0.1*95 

0.0 



—0.071? 

-0.0962 

-0.0*6* 

-0.0*05 


-0 .0330 

-0.0296 

-0.01*2 


-0.0 1*2 

-0.0296 

-O.03J0 

-0.0361 

-0.0*05 

.-J.046* 

-0.0362 

-0.0732 





0.2709 

0.4055 

0.5376 

0.660* 

0.7703 


0.9621 

3. 9*13 



n- 7765 

0. 660* 

0.5376 

0.4053 

0.2709 

0.1507 

0.8 



-ft .ftO?7 


-0.0576 

-0.0501 

-0. 0**6 

-0.0406 

-0.0333 

-0.9071 



-0.0 J33 

-0.0*06 

-0 .04*6 

-0.0501 

-0 .0376 

-0.070* 

-0.0927 

-0.0927 



0 . 1822 

ft. >7.* 

0.4059 

0.5383 

0.6613 

0. 7729 

C .973* 

0.9720 

3.927 r 

n T 07 ?r> 

0 . H75? 

0.7729 

0.6615 

0 .5383 

0.* 039 

0. 771 * 

0.1572 

0.0 




-0.0839 

-0. 0682 

-0.0590 


-0.0*73 

-0.0395 

0.007 1 


0.00T1 

-0.0395 

-0.0473 

-0.052* 

-0.0590 

-0.0682 

-0.0839 

-0.1113 





6. 9710 

ft-. 663 

0,5390 


0.7758 

0.8836 

0.9753 

3.9163 



0.7796 


0 . 3390 

0.*063 

0.2719 

0.1939 

9.0 





-0.0763 

-0.0675 

-0.0398 

-0.0333 

-0.0*22 



0.02*0 

-0.0* 22 

-0.0535 

-0.0546 

-0.0675 

-0.0783 

-0.0969 

-3 . 1 297 

-0,1297 


0^ 


0.272* 

0.4068 

0.5398 

0 .66*3 

0.779* 

0.8938 

0. 9723 

0.906* 


0.8938 

0. 779* 

0.66*3 

0* 5398 

0.4068 

0.272* 

0.1 538 

0.0 





-0.086 1 

-0.0738 


-0.0392 

-0.0*31 



0,0*2t 

-0.0*31 

-0.9592 

-0.0668 

-0.0756 

-0.0861 

-0.1096 

-6.1*73 



A.n 


0.2730 

(l.«072 

0.5*06 

0.6660 

0. 7036 

0*9057 

0.9331 

3.6903 




0.6660 

0.5*06 

0.4072 


0. 1561 





-0. 1221 

-0. 0977 

-0.0837 

-0.0733 

-0.06*3 

-0-0*18 

0.0370 

3.0330 

0.00 75 

-0.0*18 

-0 .06*5 

-0.0735 

-0 .0837 

-0.0977 

-0.1221 





0. 1 605 

0.2735 

0- *077 

0.6*1 6 

o.ee79 

0.706* 

0.916* 

0.95*6 

3.8913 




0.6679 

0.5*1 6 

0 .*077 






-0. 1 821 

-0. 1 3*4 


-0.0915 

-o.oaoo 


-0.03?* 

0.071 1 

6. 0289 

0.071, 

-0.037* 

-0. 069* 

-0.0800 

-0.091 5 

-0. 107 1 

-0.13** 

-0. 1821 

-0.1 621 




ft.??.. 

0-4062 

0.5*27 

0 .6700 

0.7939 

0.9307 

0.9*19 

3.8833 




0.6700 

0. 3*2 7 


0.27*1 

0.1632 

o.o 




-0. 1*66 

-0.116* 

-0.0991 

-0.0862 

-0.0739 

-0. 0293 

0.091 6 


0, 08 10 

-0.0293 

-0.0739 

-0.0962 

-0.0991 

-0.1164 

-0. 1466 

-0. 1969 

-0.1989 




n . > ? . r 

n.ifti? 

0 .5*38 

0.672* 

0.8002 

0.9*12 

0.9276 

3.881 3 





0.5*36 

0.4067 

0.27*7 

0. I 662 

o.o 


-0.210* 

-0.215* 


-0. 1 255 

-0. 1 066 

-0. 0922 

-0.0780 

-0, 01 79 



0.00QB 

-0.0179 

-0.0780 

-0.0922 

-0.1066 

-0,1255 

-0.1366 

-0.213* 

-0.21 3* 



n . im 

ft.?751 

0.4093 

0.3*30 

0. 67*9 

0 . 80 73 

0. 9*87 

0.91 2* 

3.0709 



0.6073 

0.67*9 

0 .3*50 

0.4093 

0.2733 

0. 1 69* 

0. 0 




-0. I 70* 

-0. 1 3*5 

-0.113% 

-0.0980 

-0.081 7 

-0.0038 

0.0931 

3.0139 

0.0R51 

-0.0036 

-0.0817 

-0.0980 

-0.1139 

-O. 1 3*» 

-0.170* 

-0.7317 



n.n 

0.1728 

0. 2759 

O.*O90 

0.3*63 

0.6776 

0.6132 

0-9327 

0.8968 

0.878* 

0. PROS 

0.9527 

0. 81 52 

0.6776 

0.3^63 

0.4098 

0.2759 

0.1728 



-o. ?• ?y 

-0.2*77 

-0. 182? 

-0.1*3* 

-0 .1 21 1 

-0. 1036 


0.81 I 7 



0.0077 

0.0117 

-0.0850 

-0. 1 036 

-0.121 1 

-0.1*3* 

-0.1822 

-0.2*77 

-0.2*77 


ft . n 


0.2766 

0.* 10* 

0.3*77 

0.6603 

0.8239 

0.9337 

0.881 0 

3.8796 


0.9532 

0.8239 

0.6805 

0.3*77 

0.410* 

0.276'' 

0.176* 

0.0 





-0.1 522 

-0.1 281 

-0.1091 

-0.0878 

0.0278 



0.0077 

0.0278 

-0.0879 

-0.1091 

-0.1 281 

-0.1322 

-6.1939 

-0.263* 

-0. 2636 


1)tA 

0.1803 

0.2772 

0.4 110 

0.3*91 

0.6836 

0.833* 

6. 930* 

0.8356 

3.8823 



O. 933* 

0.6*16 

0 .3*91 

0.411 O 

0 . 1772 

O. ISOS 



-0- 2700 

-0.2788 

-0.205* 

-0.1609 

-0.1351 

-0. 11*3 

-0.0901 

6.0*36 

0.0132 

3.00*2 

o.oooa 

0,0*36 

-0.0901 

-0.1 1*3 

-0.1 331 

-0.1609 

-0.2034 

-0. 2738 

-0.2758 


ft . rt 

0. 18** 

«.2T79 

0.4116 

0.3507 

0 .6869 

0.6*36 

0.9*4* 

0.8310 

3.083* 


0.9*** 

0.8*36 

0.6869 

0.5307 

0.4116 

0.2776 

0. 1 6«* 

0.0 


-9. ?UJ|» 

-0.2939 

-0.21 69 

-0.1696 

-0.1*19 

-0.1193 


0.0306 

0.690* 

0.00*0 

0.090* 

0.0586 

-0.0916 

-0. 1193 

-0. 1*1 9 

-0.1696 

-0.2169 

-8.2939 

-0.2939 



0.1 887 

0.278* 

0.4122 

0.332 3 

0.6903 

0.85*3 

6. 939T 

8.8976 

3. 0080 


0.9357 

0.65*3 

0.6903 

0.5523 

0.4122 

0.226* 

0.1887 

0.0 




-0.2263 

-0.1 761 

-0.1 *63 

-0. 12*0 

-0.0922 

6.6723 



0. 0815 

0.0723 

-0.092* 

-0.1 2*0 

-0.1*65 

-0.1781 

-0.2283 

-0. 3088 

-0. 3066 


0.0 

0.1933 

0.2790 

0.* 128 

0.35*0 

O .69*0 

0.865* 

0.92*« 

0.0237 

3.8829 




0.6940 

0.53*0 


0.2790 

6.1933 

0.0 


-0. 123* 

-0.323* 

-0.2395 

-0. 1865 

-0.1331 

-0.1283 

-6.0919 

6.08*2 

6. 0768 

0.8086 

0.07*8 

0.08*2 

-0.3919 

—9 . | 285 

-0.1331 

-0. 1 863 

-0.2399 

-6. 323* 




Figure A-3. Continued. 
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SIDCVA9H INFORMATION FOM POSITIVE T-AXlS AT T»« * STATION OF 0*6000 3ENI-3PAN3 


0.1 766 
O. I *6 6 
O * 07 6 A 
0 . 0*66 
-0.0232 
-0.0 7 32 
-0.1232 
-0,1 732 
-0*2232 
-0.2732 
-0.1232 




■ 0.4031 


0*2036 
0.1336 
0*1036 
0*0336 
0*0036 
- O* 0464 
-0.0464 
-O. 1464 
-0. 1464 
-0*2464 
-0.2404 


-1.36 

-0*40 

0.03 


■ 6A 
.04 

► 03 

► 77 


V • 0*3371 
t A MG LI 

0*2043 -1*4! 

0*1343 -1.3! 

0*1043 -1.33 

0.0343 -1.34 

0.0043 -0.81 

-0*0403 0,01 

-0.0403 0.4: 

-0.1403 1*4! 

-0.1403 1.31 

-0.2403 1.31 

-0.2403 1.«« 



V - 0.7606 
2 ANGLE 

0.2170 -0.96 

0.167Q -0.61 

0.1170 -0.64 

0 . 0670 —9.36 

6.0176 -0.26 

-0.0330 0.23 

-0.0830 0.62 

-0.1330 0.72 

-0.1330 0.64 

-0.2330 0.61 

-0.2030 ' 0.34 


004NV4SH INFORMATION F 0 R POSITIVE T-AXIS AT THE X STATION OF 0.6000 SEMI-SPANS 



-0.1464 2.41 
-0.2464 2. 28 
-0.2464 2. 13 


T » 0.3371 
z ANGLE 

0.2043 3.06 

0.1343 3.21 

0.1093 3. 3 A 

O.0S43 3.43 

0.0093 3.43 

-0.0403 3.33 

-0.0405 3.14 

-0.1403 3.06 

-0.1403 2.46 

-0.240 3 2.83 

-0.2903 2.64 




-0.2333 4.01 

-0.2830 3.61 



Figure A-3. Continued. 
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APPENDIX B-WAKE PROGRAM 


User Instructions 


The program is written in FORTRAN IV and is designed to run in double 
precision on an IBM 370-165 computer* Using a finite difference technique 
this program evaluates the velocity profile downstream of a wing trailing 
edge given the upper and lower surface profile shapes at the wing trailing 
edge. All the input information can be obtained from the two-dimensional 
airfoil program in Reference 4. Using 1601 points in the n or normal direc- 
tion with an initial stepsize of .00005 in £(x/c) anc | an increment in step- 
size of .05, an average execution time of 2 minutes and 42 seconds was 
required to integrate the boundary layer equations downstream two wing chord 
lengths. The program requires the specification of the following input data: 


Card Variable 


Number 

Name 

Variable Description 

1 

TITLE 

The 80 character array which is used as a header to identify 
output. Termination of execution is achieved by following 
the last set of profile data to be analyzed by a title card 
having only the word END in the first three spaces. 

2 

DETA 

Stepsize in the n-di rection . A value of .00125 was found 
satisfactory for most cases, and using this value with a 
constant stepsize, the total number of n-points, NPTS, can 
be estimated using the equations given in the text. 

2 

DX 

Stepsize in the X-di rection. Best results were obtained 
for cases analyzed in this report using a very small step- 
size, .00005 (specified in chords), and unequal spacing 
using the variable DXSTEP on card 3. 

3 

X(1 ) 

The station in chords at which the initial velocity profile 
solution is assumed. A value of .01 was found satisfactory 
for most wake solutions. 

3 

CEDDY 

A constant used to evaluate the eddy viscosity model 
required to solve for the velocity profile. A value of 
.03 correlated best with experimental data. 

3 

T 

The thickness at the spanwise station of interest in chords. 

3 

CD 

The drag coefficient at the spanwise station of interest. 

3 

CDMIN 

The minimum drag coefficient of the wing section at the 


spanwise station of interest. 


Ill 



Card Variab le 

Number Name Variab le Descri ption 

3 DXSTEP Steps i ze increment in the X-di recti on. The steps ize in the 

X-direction is increased by the factor DX , ,, = 

' step n+1 

(DX s _|_ e p n )(1 + DXSTEP). A table is given in Figure 30 which 

will help in determining the X-va I ue at any station using 
unequal step increments of .05, .075, and .10. 

4 NXPRNT The NXPRT station numbers at which the user desires to have 

the velocity profile printed. The maximum number of stations 
which may be printed is 20 unless the solution is printed at 
every station (NXPRT=0 and all values of NXPRNT are specified 
as 0). The station number corresponds to a specific value 
of X downstream; if equal stepsize in the X-direction is 
used then the NXPRNT values at the desired X-station are 
easily calculated. However, for unequal stepsize which is 
the more useful case Figure 30 will help in determining the 
step number nearest the X-station at which the velocity 
profile is desired. The NXPRT values of NXPRNT are read 
using a 2014 format. If only two stations were of interest 
then only columns 1-4 and 5-8 would contain integers 
representing the station numbers. The stations must be 
specified in the order of increasing station number. 

5 DELTAU The boundary layer thickness in chords at the trailing edge 

on the upper surface at the station of interest. 

5 DELTAL The boundary layer thickness in chords at the trailing edge 

on the lower surface at the station of interest. 

5 HU Form factor on the upper surface at the station of interest 

= boundary layer displacement thickness divided by the 
momentum thickness at the trailing edge. 

5 HL Form factor on the lower surface at the station of interest 

= boundary layer displacement thickness divided by the 
momentum thickness at the trailing edge. 

Cards 1 through 5 represent a complete data set for the WAKE Program. 
Several data sets may be executed for a single program compile by placing 
subsequent data sets behind the first. The entire program execution is 
terminated when the END title card is encountered. The format specification 
for the above data is given in Figure B-1. A sample data set Is given in 
Figure B-2, and the output of this data set is presented in Figure B-3. 
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Figure B-1 . Format specification of input data for the WAKE program 














Figure B-2. Example data set for the WAKE program 
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Program Listing 


C444 , USING A FINITE DIFFERENCE METHOD THIS PROGRAM CALCULATES THE 2-0 

C444 'velocity profile In The wake of an airfoil given the boundary 
C**« layer information on the upper and lower surfaces of the airfoil 
C444 AT THE TRAILING EDGE* 

IMPLICIT REALF0I A-H.O-Z) 

COMMON A.S.C.ft. X ( 300 ).F(200l.3).W(200l.3)« WOLD ( 20011 .CAPE< 2001 ) . SM 
1 ALEC 2001 ) *E TAT 2001 ) • WALBC » OEBC • AK. OX l I .DXI2.DX13.BB.C1 ■ I TER • NX S* n , 
INPTStNEND 

DIMENSION T|TLE(20) .NXPRNT ( 20 ) 

DATA EN0/AHEN0 / 

IREAD-I 
!PUNCH»2 
t MR I TE = 3 

C*44 READ TITLE CARO 

5 RE AO TIREAD.10) ( T I TLE( I I • I* 1 • 20 * 

10 FORMAT (20A4> 

IF (TITLE! I)«E0. END) GO TO 170 
C44 4 INPUT BASIC CONSTANTS 
C*»« NPTS*NUMBER OF ETA POINTS OESIRED 

C«** NET APT "NUMBER OF STEPS BETMEEN ETA POINTS PRINTED 
CM* NX STEP* NUMB EM OF STEPS IN x DIRECTION 

c*** NXPRT-NUMBER OF STATIONS AT WHICH SOLUTION IS To BE PRINTED 
C444 CNXPRT»0 PRINTS EVERY SOLUTION AT EVERY STATION) 

C444 (THERE IS A MAXIMUM OF 20 STATIONS) 

C444 (THE FIRST STATION IS ALWAYS PRINTED) 

CM» (THE STATIONS MUST BE SPECIFIED IN ORDER OF INCREASING 

C«*« STATION NUMBER) 

C • • • T*TH ICKNESS TO CHORO RATIO OF MING 

C444 912=WAKE WIDTH FACTOR 

C444 UMUDEL*VEL0CITY OEFICIT FACTOR 

CAM* CK*CONSTANT USED IN EDDY VISCOSITY MODEL I 

C444 A K* PROPORT I ONAL I TV CONSTANT FOR UNEVEN STEP SIZE IN ETA DIRECTION 
CAM* EPS* ACCURACY FOR SOLUTION CONVERGENCE 
C44* OETA*INITIAL STEPSIZE IN ETA QfRECT ION 

C444 0X*STEP5IZE IN X DIRECTION (INPUT AS RATIO TO CHORD LENGTH) 

READ ( I READ* 15 I NPTS. NETAPT • NX STEP , NXPRT . AM .EPS.DE TA . OX 
IS FORMAT (4U0«4F|0.5) 

READ (IREAD.20) XI 1 ) . CEOOY « T .CO.COM IN* DXSTEP 
20 FORMAT ( SF 10*5) 

NXPRT 1*NXPRT 

IF ( NXPRT.EO.O)NXPRTl=t 

READ (IREAD«2S) (NXPRNT! I) .1=1 •NXPRTlI 
25 FORMAT (2014) 

READ ( I READ* 20 I OELT AU. DELT AL.HU.HL 
NEND-NPTS-l 
I NXPRT*! 

C444 SET OUTER EDGE BOUNDARY CONDITION AT 1.0 
OEBC* ! *0D0 

PU*2.D0/( HU-1 .00 I 

PL *2.00/1 ML- 1*00) 

612*. 500 • t • 7900 4 4PU*0EL TAU*. 7 5D0 **PL 40EL T AL I 
UHUOCL* I CO-COM IN ) /CO 

Cl"T«CD«tOO.OO/(0l2*CEOOY*( 1 .DQ-UMUDEL) ) 

NXS*l 

NMIO*(KPTS-l 1/241 

WRITE (| WRITE. 30) (TITLE Cl) .1*1 .20) . NPTS .NMID. NETAPT • NXSTEP. NXPRT t 
I (NXPRNT! I )• 1*1 .NXPRT I I 

30 FORMAT ( IHl • //Sx .20A4 • //5X.9HNPTS * . I4./SX .9HNMID = . 14./5X.9 

1HNETAPT ■ .I4./5X.9HNXSTEP - . I 4 ./5X , 9HNXPRT * . I 4 ./5X .9MNXPRNT = 

1 .2014//) 


WRITE ( I WRITE. 35) AK . EPS. DE TA. DX , X ( 1 ) .CEDDY. B l 2. T. UMUOEL .CO .OXSTEP 
1.C1 

■ .DIS.8./5K. 9HEPS * .D15.B./5X.9HDETA * 

* .0|5*8./5X.9HX( 1 ) = .015.8./5X.9HCE00Y ■ 

= .D15.8./5X.9HT = .01 5 .8, /5X .9HUMU0EL * 

= .015.S./5X.9H0XSTEP * , 01 5. 8. /5X.9HC 1 * 


3S FORMAT (5X.9HAK 

1 .015. 8. /5X.9HOX 
1 .015.6. /5 X • 9HB 1 2 

1 . 015. 8. /5X.9HCD 
I.D1S.6//) 

WRITE (twRITE.40) DELT AU .DEL TAL .HU. HL 


•015.8. /5X.9HHU 


40 FORMAT (5X.9HDELTAU = . 01 5 . 8./5X.9HDELT AL 
1 .015.6. /SX.9HHL * .015.6//) 

C*44 CALCULATE E T A ( ( ) 

ETA(NMID) *0 • DO 
ETAINMI0-1 >=-DETA 
ETA! NM 10+ 1 l*0ET A 
J * NM1 D— 1 
JJ=NMI D+2 
DO 45 I=JJ, NPTS 
J*J-1 

ETA( I ) *£TA( 1-1 ) 4 AK ♦ ( E T A( I-|)-ETA( 1-2)1 
45 ETAl J)=-ETA( I ) 

DO 50 1=1. NPTS 
F ( 1.21=0.000 
F< I .31*0.000 
W(| . 2) *0 . 000 
50 W( I * 3)*0.000 
OO 55 1*2*300 
IF ( 1 • EO • 2 ) GO TO 55 
OX=DX*< I .D04DXSTEP) 

55 X ( I |«X( 1-1 )*DX 
c*4* INITIAL PROFILE GUESS 
CU=( HU-1 .00) /2.D0 
CL=(HL-1 .00 )/2 .DO 

OELTAU=OELTAU/OSO«T(T*X(NXS)41 00. 004 CO) 
OELTAL*OELTAL/DSQ«TtT4X(NXS> 4t00«00 4CD I 
CALL PROFIL(OELTAU.OELTAL.CU.CL.OETA.NMID) 

00 60 1*1 .NPTS 
60 WOLDC I > = W( I. 1 I 
65 IF (NXS-2) 75.75.70 

70 DX l 1* ( ( 1.000/! X( NXSt-X! NXS-ll) ) 4 ( 1 .000/! X( NXS) -X( NXS-2 I ) I ) 

DX12=( X(NXS)-X (NXS-2) ) /( ( X( NXS )—X( NXS- I ) ) *( X(NXS-1)-X( NXS-2) I ) 
0X1 3=( X! NXS ) -X ( NXS— 1 ) )/(( X( NXS I -Xl NXS-2 ) ) 4( X( NXS-1 )-X( NXS-2))) 
75 CONTINUE 

F ( NM 10 . 1 )*0 .DO 

C 4 * 4 START SOLUTION AT FIRST STATION 
ITER- I 
80 N* 1 

CALL GETFINMIDI 
IF (NXS.EQ. 1 ) GO TO 110 
as n=nm 

C • • * CALCULATE MATRIX COEFFICIENTS 
CALL A0CR 
CALL COEFF 

IF (N-NENDI 85.90.90 
C444 CALCULATE NEW W ARRAY PROFILES 
90 CALL SOLVE 

C444 TEST FOR CONVERGENCE OF SCCESSIVE PROFILES 
OO 95 1=2 . NEND 

[F (OABSI 1.000-1 MOLD! I )/W( 1. t) ) 1-EPSl 95.95.100 
95 CONTINUE 

C * 4 4 CONVERGENCE ATTAINED 


FIN 61 
FIN 62 
FIN 63 
F(N 64 
FIN 65 
FIN 66 
FIN 67 
FIN 68 
FIN 69 
FIN 70 
FIN 71 
FIN 72 
FIN 73 
FIN 74 
FIN 75 
FIN 76 
FIN 77 
FIN 78 
FIN 79 
FIN 80 
F|N 81 
FIN 92 
FIN 83 
FIN 84 
FIN 85 
FIN 86 
FIN 87 
FIN 88 
FIN 89 
FIN 90 
FIN 91 
FIN 92 
FIN 93 
FIN 94 

FIN 95 
FIN 96 
FIN 97 
FIN 98 
FIN 99 
FIN 100 


FIN 103 
FIN 104 


FIN 107 
FIN 108 


FIN 115 
FIN 116 



116 


<«o to no 

C444 CONVERGENCE NOT ATT A INEO- I TERAT I ON SOLUTION 
100 CONTINUE 

IF I 1 TER sGT *99 I GO TO 160 

(TER-VTER+I 

DO 105 1*1, NPTS 

WOLD| I l>l( I,|) 

105 CONTINUE 
GO TO 00 

CM* PRINT CONVERGED SOLUTION 
no if (NXPRT.eo.oi go to iis 

IF (NXS.EO.l ) GO TO 115 

IF | NX S, EQ ■ NXPRNTI 1 NfcPRT ) ) GO TO 115 

GO TO 150 

IIS WRITE ( I WR( TE, 1351 ITER 

WRITE ( I MR I TE, 1 20) NXS.X(NXS) 

120 FORMAT I/I0X.23HRESULTS FOR X/C STAT (ON* I 4,5X,4HX * ,015.0) 

WRITE (IWRITE.125) C T I TLE< I ) , 1 *1 ,20 ) 

125 FORMAT ( /• I X,20A4./> 

WRITE ( IWRITE, I 40 > 

C YET AxOSORTl TAXI NXS >41 00 ,00* CO ) 

l-npts+netapt 

00 130 1*1 ,NPTS«NETAPT 

L*L-NETAPT 

YETA*ETA«L)*CVETA 

WS0*W(L,l l»W<L.l > 

WRITE 1IWRITE.IA5) ET A ( L I.W(L.l) *F I L * t ) , VETA • WSQ 
130 CONTINUE 

IF ( NXS* N£, 1 ) I NXPRT*INXPRT 4 1 

|F < NXPRT,NE,0*ANO. INXPRT,GT.NKPRT) GO TO 5 
135 FORMAT I 1 HI • ✓ , I OX, 26HC0NVERGENCE ATTAINED AFTER, 14# 1 IH ITERATIONS) 
I 

140 FORMAT </6X«3HETA, I2X«6HU/UDEL*13X • t HF . I 4X ,3H*/C ,9X , 1 i HI U/UDEL ) • 4 2 
I I 

145 FORMAT I 5016,0 ) 

C444 SET UP FOR NEXT X STATION 
ISO 00 155 1*1. NPTS 
W (1 ,3>*W<I.2> 

F ( t*3)*F( |,2) 

Wf I ,2 )*Wf I,| ) 

Fit *2) *F< 1,1 ) 

WOLDI 1 l-wt t, | ) 

155 CONTINUE 
NXS-NX541 

|F <NXS-Z,GE,NX5TEP ) GO TO 5 
GO TO 65 

160 WRITE I I MR I TE , 1 65) NXS 

165 FORMAT </,5X»30MlTER EXCEEDS T5 AT STEP NXS * .141 
GO TO 5 
170 CALL EXIT 


FIN 121 
FIN 122 

FIN 123 

FIN 124 
FIN 125 
F|N 126 
FIN 127 
FIN 120 
FIN 129 

FJN 130 
FIN 131 
FIN 132 
FIN 133 
FIN 134 
FIN 135 
FIN 136 
FIN 137 
FIN 130 
FIN 139 
FIN 140 
FIN 141 
FIN 142 
FIN 143 
FIN 144 
FIN 145 
FIN 146 
FIN 147 
FIN 140 
FIN 149 
FIN 150 
FIN 151 
FIN 152 
FIN 153 
FIN 154 
FIN 155 
FIN 156 
FIN 157 
F|N 150 
FIN 159 
FIN 160 
FIN 161 
FIN 162 
FIN 163 
FIN 164 
FIN 165 
FIN 166 
FIN 167 
FIN 160 
FIN 169 
FIN 170 


END 


FIN 171 


SUBROUTINE GETF(NMIO) GET 1 

C444 THIS SUBROUTINE INTEGRATES THE DERIVATIVE OF F TO OBTAIN F GET 2 

IMPLICIT RCALWBf A-H.O-ZI GET 3 

COMMON A*B,C.R,X<300),F<2001, J), W< 200 1 . 3 ) , WQLDI 200 1 ) , CAPE1 200 I I.SM GET 4 

I ALEC 2001) .ETAI 2001 ) , tf ALBC ,OEBC • AX.DX 11,0x12,0X1 3* BB, Cl » I TER* NX $« N • GET 5 

1M**T S* NEND GET 6 


FINL*0,D0 GET 7 

FJNU*0,00 8 

JJsNNtDM GET 9 

J-NMID GET 10 

JJJ-J-2 G ^T II 

NPTSM\*NPTS-t GET t2 

DO 5 I=JJ,NPTSM1 GET 13 

F(NU=F1NU*( WOLDI t)+tfOLD< I-l))*tETAI 1 ) -ET A I I - 1 1) /2. DO GET 14 

5 F< I, | )-FlNU GET 15 

DO 10 IsJ.JJJ GET 16 

jwj-I GET 17 

F INL«FI*4.4|WOLO(J)4WOLOU4l) I4IETAC J)-ETAI J41 D/2.D0 GET 10 

10 F|J.I)*FINL GET 19 

FI NPTS, 1 ) =F ( NENO , 1 )*■< WINPTS.1 >4W(NEND.l )) *| ET A I NPTS 1 -E T A( NPTS- l ) )/ GET 20 

12,00 GET 21 

FI |»l)aFl2,l)+|M(l«l) + w(2,l))* IET A( ! )— ET Al 2) ) /2.D0 GET 22 

RETURN GET 23 

END GET 24 


SUBROUTINE A0CR , ABC 

C444 THIS SUBROUTINE COMPUTES A,B,C, AND P MATRIX COEFFICIENTS ABC 

IMPLICIT REAL*0l A-H.O-2) ABC 

COMMON A,B,C.R.X ( 300 ). FI 2001 • 3 ) • V ( 20 0 I • 3 ) . WOL D I 200 1) . C APE I 2001 1 • SM ABC 
1 ALE I 2001 ) ,E TaI 2001 ) , WALBC • OEBC. AK, Dx 1 1 ,DX! 2, DX1 3 , SB, Cl « ITER • NX S, N , ABC 
1 NPTS. NEND ABC 

O t MENS ION ALPHAI2) ABC 

Dlx|ETA(N4| ) -ETA|N))4AK*AK*lETA|N)-ETA<N-ll ) ABC 

02-1 ETAI NM >— ETA(N) ) 4424AK •{ ETAC N )-ET Al N- I ) ) 442 ABC 

C*** COMPUTE ALPHAS FOR MOMENTUM EOUATION ABC 

IF (NXS-2) 5*10,15 ABC 

5 ALPHA 1 1 )*CI *F <N, I )/2 ,D0 ABC 

ALPHA! 2)*0. DO ABC 

A* f 2«0D0/D2 + ALPHA! 1 ) /D 1 1 ABC 

B=-Z .0004(1 ,0D0*AK)/O2-( 1 . 0D0— AK4AK ) 4 ALPHA I 1 ) /O 1 ABC 

C=2.0D0*AK/O2-ALPMA( I ) 4AK4AK/01 ABC 

R*0 .DO ABC 

GO TO 20 ABC 

10 ALPHA! l)x(F|N, 1)/2«D0+X(NXS) 4|F(N.l )-F(N,2) )/( XI NXS I -X < NXS-1 ) ) )4Cl ABC 

ALPHA! 2I»-X(NXS14W0L0(N)4CI ABC 

A = ( 2.000/D2 *■ ALPHA! 1 1/DI ) ABC 

8 --2 .0D04U .0O0*-AK)/D2-| l,OD0' AK4AK) 4 ALPHA 11 )/Dl*-ALPHA 12) /(X(NXS)^ ABC 
IK(NXS't)) ABC 

C“2, 0DO4 A K/D2- ALPHA I 1 )4AK4AK/D| ABC 

R* ALPHA I 2)*W( N,2)/| Xl NXS) -X< NX S- 1 I ) ABC 

GO TO 20 ABC 

15 ALPHA! 1 ) =f F( N, 1 )/2 ,004X| NXS) * (Ox I 14F(N,1 I — DXl 24F< N, 2 ) *DXl34Ff N ,3) ) A0C 

J>*Cl ABC 

ALPHA! 2) — X( NXS I4WOLOIN) 4CI ABC 

Ax ( 2, 0D0 /0 2-4 ALPHA! I ) /D1 I ABC 

0*1-2,0004! 1 ,0D04AK.)/D2-n ,0O0-AK4AK)4ALPHA« I ) /0 1 4 ALPH At 21 40X 1 1 ) ABC 

C-2. 0004 AK/02- ALPHA! 1 )*AK*A*/Dt ABC 

R -ALPHA I 2)40X1 24 «| N ,2 ) -ALPHA I 2 ) 4DX 1 34W ( N , 3) ABC 

20 RETURN ABC 

END ABC 


to 
I 1 
12 
13 
1 4 

15 

16 
17 
10 

19 

20 
21 
22 

23 

24 
29 
26 
27 
20 

29 

30 

31 

32 

33 

34 

35 


SUBROUTINE COEFF 


COF 


! 


i 

i 


1 17 



C44* COEFF COMPUTES RECURSION COEFFICIENTS FOR TRIDIAGIONAL COF 2 

C444 MATRIX SYSTEMS COF 3 

IMPLICIT REAL*a< A-H.O-Zl COF 4 

COMMON A.B.C.R. X( 300) .F( 2001 *3) .M( 2001 .3) «WOLD(200l) «CAPE< 2001 >.SM COF 5 

I ALE 1200 1 ) »ET At 2001 I .MALBC. OEBC. AK. DX I I .DX 1 2.DX » 3.0B.C I * I TER • NX S. N. COF 6 

INPTS.NENO COF 7 

IF IN-2) 5.5.10 COF a 

5 CAPEIN1*— A/0 COF 9 

SMALE<N)*(R-C«OEBC)/a COF 10 

GO TO IS COF II 

10 C APE I N )*- A/ 1 0*C*CAPE( N— I I) COF 12 

SMALE<N) = (R-C*SMALE(N-n I/* a+C*CAPE(N-m COF 13 

15 RETURN COF 14 

ENO COF 15 


SUBROUTINE SCLVE SOL I 

C 4 4 4 THIS SUBROUTINE SOLVES A TPIDIAGIONAL MATRIX SYSTEM BY SOL 2 

C444 SACK-SUBSTITUTION SOL 3 

IMPLICIT REAL48I A— H.O-*Z> SOL 4 

COMMON A.B.C«R»X(300) *F < 2001 .31 « V( 2 00 1 .31 .VOL 01 2001 1 ♦ CAPE ( 2001 ) * SM SOL 5 

1 ALE (2001 > • ET AC 200 1 ) .MALBC. OE SC. AK* DX 1 l.DK 12.0X1 3. 8B.C 1 • I TER. NAS. N. SOL 6 

INPTS.NENO SOL 7 

N*NM SOL 8 

■|N.l)=OEBC SOL 9 

5 N*N-1 SOL 10 

M(N. n*CAPE(N)4M<N+l . 1)4SMALE( N) SOL 11 

IF ( N— 2 ) 10.10.5 SOL 12 

10 Ml 1 .1 >*OE0C SOL 1 3 

RETURN SOL 14 

END SOL 15 


SUBROUTINE PRCF IL C DELT AU. OELT AL. CU. CL* BETA , NM ID ) PRO 

C • * 4 THIS SUBROUTINE EVALUATES THE INITIAL VELOCITY PROFILE FOR THE PRO 

C44* FINITE DIFFERENCE SOLUTION GIVEN THE BCUNOAOY LAYER THICKNESS AND PRO 

C •* 4 FORM FACTOR FOR BOTH UPPER ANO LOMER SURFACES PRO 

IMPLICIT REAL48I A-H.O-21 PRO 

COMMON A.B.C.R.Xi 300 >.F(20Ol .3I.MI 200 1 , 3 ) ,MOLD( 200 1 ) .CAPE! 2001 ) .SM PRO 
I ALE (2001 ) .ETA I 200 1 1 • MALBC. OE0C • AK. DX 1 1 . DX 1 2. DX 1 3 . BB. C 1 . 1 TER *NA S. N . PRO 
INPTS.NENO PRO 

NSTARTsNMlD PRO 

DELETA«-D€TA PRO 

00 IQ I*N5T ART .NPT5 PRO 

DELETA-DELETA4DETA PRO 

IF (OELETA.'GT.DELTAUI GO TO 5 PRO 

MM. I l*(OELETA/OELTAU) 44CU PRO 

GO TO 10 PRO 

5 M ( I. lt*l.D0 PRO 

to CONTINUE PRO 

deleta=o.oo pro 

NNENO»NSTART-l PRO 

NCNT»NST ART PRO 

00 20 1*1. NNEND PRO 

NCNT=NCNT-1 PRO 

OELETA*OELETA*OETA PRO 

IF IOELETA.GT.OELTAL) GO TO 15 PRO 

MINCNT .1 )*IDELETA/DELTAL)44CL PRO 


10 
i t 
12 
1 3 
1 4 

15 

16 
1 7 
18 

19 

20 
21 
22 

23 

24 

25 


I 


GO TO 20 PRO 26 

15 M( NCNT .11*1 .DO PRO 27 

20 CONTINUE PRO 28 

RETURN PRO 2R 

END PRO 30 





Somple Output 


I 


/ 


JOUlCCaSKI AIfiFOlL/A = 0/HN«. 21/TREE TRANSIT rON/Ma.OS/CEODY*.OV 


W*TS 

NMIO 

NETAPT 

nxstcp 

nxopt 

NAPRNT 

EPS 
DETA 
OX 
X< 1 ^ 

ceooy 


teoi 
no i 
20 

ISA 

113 120 136 157 
0*100000000 01 
O. 1 00000000—04 
0. 12500000D-02 
0.500000000-04 
O. I OCOOOOOD-O I 
0. 300000000-0 1 


PI 2 
T 

U MODEL 
CO 

OXSTEP 

Cl 


0.604QA6 I 40-02 
0. 116000000 00 
0.0 

0. I 03780000-01 
O.SCOOOOOOD-OJ 
0.674662270 03 


OELT AU 
OELTAL 

HL 


0.24267000D-01 
0.242070000-0 1 
0.141390000 01 
0. 141390000 01 


JOUKCtfSK I 


CONVERGENCE ATTAINED AFTER 1 
RESULTS FOR K/C STATION I 
AIPF0IL/A>0/BN».21/PPE£ TRANSI 


I TERAT IONS 
X * 0.100000000-01 

TI0N/Ma.05/CCO0V=.03/ 


0.100000000 01 
0.975000000 00 
0.950000000 00 
0 .925 OOOOO 0 00 
0.90000 OO OD 00 
0.0 75000000 00 
0.660000000 00 
O • 62 50000 OD OO 
0 . 0000000 OD 00 
0.775000000 00 

0.750000000 OO 
0.725000000 00 
0.700000000 00 
0.675000000 00 
0.650000000 00 
0.625000000 00 
0.600000000 00 
0.575 0 OOOOD 00 
0.550000000 00 
0.525000000 00 
0.500000000 00 
0 . 4 750*1 GOOD 00 
0.450000000 00 
0.425000000 00 
0.4000000GD 00 
0. 375000000 00 
0.350600000 00 
0.325000000 00 
0.300000000 00 
0.275000000 00 
O.25OOQ0OOD 00 
0.225000000 00 
0.200000000 00 
0.175000000 00 
0.150000000 00 
0.125000000 00 
0.100000000 00 
0. 750000000-0 1 
0. 5 OOO OOOC 0-0 l 
O .250000000-0 1 
0.0 


— 0.25000 0000—0 I 
-0.500000000-01 
-0. 750000000-0.1 
- 0.100000000 00 

-0. I 25000000 00 
-0.150000000 00 
-O.I750QOOOD OO 
-0.200000000 00 
-0.225000000 00 
-0.250000000 00 
-0.270000000 00 
-O.JOOOOOOOD 00 
-0.325000000 00 
-0.350000000 00 
-0.375000000 00 
-0.400000000 00 
-O.425OOO0OD 00 
-0.450000000 00 
-0.475000000 00 
-0.500000000 00 
— 0. 5250 OOOOD 00 
-0.550000000 00 
-0.575000000 00 
- 0.600000000 00 

— 0 .625000000 00 
—0.65000 0000 00 

— 0 . 6750Q 00 Oo 00 
-0.700000000 00 
-0.725000000 00 
-0.750000000 00 
-0.775000000 00 
-C.SOOOOOOOO 00 
-0.925000000 00 
-0.650000000 00 
-0.875000000 00 
-0.900000000 00 
-0.925000000 00 
-0.950000000 00 
-0.9T500000D 00 
- 0.100000000 01 


U/UOEL 

0.100000000 01 

O.IOCOOOOOO 01 

0.100000000 01 
0. I 000 0 0 0 OD 01 
0.100000000 01 

O. 10000000D 0 1 

o.ioooooooo 01 

O.1OOOOQO0D 01 

0.100000000 01 
O.IOOOOOOOO 01 
O.IOOOOOOOO 01 

o . 100000000 01 
0.100000000 01 

0.994263640 00 
C. 986528530 00 
0.976553550 00 


0 *961 61 2660 00 
0.95300524D 00 
0.943874390 00 
0.934391910 00 
0.924525670 00 
0.9142366 IO 00 
0.903497850 00 
0 . 992 2232 70 OO 
0.960365740 00 
0.867904840 00 
0.854695640 00 
0.94065437D 00 
0.825652180 00 
0.809526210 00 
0.792066060 00 
0.772992690 00 
0.751924010 00 
C. 726315100 00 
0.70I34670D 00 
0.66969527D 00 
0.6 30968080 00 
0.580201600 00 


O . 66 C 691 600 00 
0.955891900 00 
0.630981800 00 
0.905881900 00 
0.780881600 00 
0 . 75566 ! BOD 00 
0.730861800 00 
0.705981800 00 
0.68098160D 00 
0.655681900 00 
0.630691800 00 
0.605881800 00 
0*580981900 00 
0.5SS9362ID 00 
0.531 I 7 56 3D 00 
O .5066 I I 780 00 
0.482250290 00 

0.456098020 00 
0.434162150 00 
0.410450460 00 
0.386971360 00 
0 » 36 37 34 0 6D 00 
0.34074059D 00 
0.316026000 00 
0*295570490 00 
0.273419610 00 
0.251564560 00 
0.230030440 CO 
0.204936710 00 
0. 188005730 00 
0.167S6347D 00 
0.147540530 00 
C. 127973600 00 
0.109907500 00 
0.903985190-01 
O. 7251 9646D-01 
0 .55 J6 9994D-0I 
0 . 390931 270-01 
0. 239 1941 40-01 


T/C 

0.34994J420-01 
0.341 194840-01 
0.332446250-01 
0. 323697670-01 
O. 31 4949CAO-OI 
O. 306200500-01 
O. 297451 910-01 
0.289703320-0 l 
0.279954 740-01 
0.271 206150-01 
0.262457570-01 
0.253708980-01 
0.244960400-0 I 
0.23621 1610-01 
0.227463230-01 
0.2187 14640-01 

0. 2 09966050-0 1 

0. 201217470-01 
0. 192466880-01 
0. 163720300-01 
0. 1 74971 710-01 
0. 166223 130-0 I 
0. I 57474540-01 
C .149725960-0 I 
0. 1 39977370-01 
0. 131 228760-01 
0 . 1 22480 20 0-01 
0.1 137316 10-01 
0.1 04983G30-0! 
0.962344420-02 
0.8 74658560—02 
0. 76737270D-02 
0 .699686050-02 
0.6 1240099D-02 
0.524915140-02 
0.437429260-02 
0.34994 34 2D — 02 
0 .262457570-02 
0. I 74971710-02 
0 .674658560-0 1 


C. 0 

0.50266762D 00 
0.560201800 00 
0.630968060 00 
O . 6696 952 7D 00 
0.701 346700 00 
0.728315100 00 
0.751924010 00 
0.7 72 992690 00 
0 .792066060 00 
0.609526210 00 
0.625632160 OO 
0.840654370 00 
C. 654695640 00 
0.66 79 04 64 D 00 
0.980385740 00 
0.89222J270 00 
0.90348765D 00 
0.914236610 00 
0.92452S670 00 
0.934J91910 00 
0.943B7438D 00 
0.953005240 00 
0.961812660 00 
0.970321450 00 
0.976SSJ55D 00 
0.966520530 00 
0.994263840 00 

o. 1 00000000 01 

0. 100000000 0 I 

0. 1 00000000 01 

O.IOOOOOOOO 01 
O.IOOOOOOOO 01 
o.ioooooooo 01 
0.100000000 01 
0.100000000 01 

O.IOOOOOOOO 01 

0 . 1 oooooooo 0 1 

O.IOOOOOOOO 01 

o.ioooooooo 01 

0.100000000 01 


-0. 102957160-01 
-0 . 239 1941 40-01 
-0. 390931 270-01 
-0.553699940-01 
-0.725196480-01 
-0.903965190-01 

— 0 • 1 0 090 7 SO O 00 
-0.127973600 00 
-0. 147540530 00 
-0.167563470 00 
-0.166005730 00 
-0.206836710 00 
-0.230030440 00 
-0.251564560 00 
-0.273419610 00 
-0.295576480 00 
-0. 31602600D 00 
-0.340748590 00 
-0.363734060 00 
-0.386971360 00 
-0.410450460 00 
-0.434162150 00 
— 0 • 456096020 00 
-0.462250290 OO 
-0.506611760 00 

— 0.531 I 758 3D 00 
-0.555936210 00 
-0.560661800 00 
-0.605881800 00 
-0.630881800 00 
-0.655881800 00 
-0.680861800 00 
-0.705881800 00 
-0.730881800 00 
-0.755881800 00 
-0.760861800 00 
-0.805881800 00 
—0.630681800 00 
-0.655881800 00 
-0.880881800 OO 


0.0 

-0.874656560-03 
-0. 1 7497 t 710-02 
-0.262457570-02 
-O. 34994 3*20-02 
-0.4 3 7*29280-02 
-0.52*91 514D-02 
-0.612400990-02 
-0.699886 850-02* 
-O. 7 67372700-02 
-0.874856560-02 
-0.962344420-02 
-0.1 04983030-01 
-O.ll 37 31 610-0 I 
-0.122460200-01 
-0. 131228780-01 
-O . 1 39977370-01 
-0. 140725960-01 
-0.157474540-0! 
-0.166223130-01 
-O. I 74971 710-01 
-0.103720300-01 
-O. 192468680-01 
-0. 201217470-01 
-0.209966050-01 
-0.2 18714640-01 
—0.22 7463 ? JO— 0 1 
-0.23621 I 610-01 
-0.244960400-01 
-0.253708980-01 
-0.262457570-01 
-0.271 206150-01 
-0.2 79954 740-01 
-0.288703320-01 
-0.297451 910-01 
-0.306200500-01 
-0.314949090-01 
-0.323697670-01 
-0. 332446250-01 
-0.341 194840-01 
-0. 349943420-01 


IU/UDEL l**2 
0.100000000 01 
3.1 OOOOOOOO 01 
O.IOOOOOOOO 01 

0. 100000000 01 
0.100000000 01 

". IOOOOOOOO Qt 

o.ioooooooo 01 

0 . I 0 0000 OOD 01 
0. IOOOOOOOO 01 
O . IOOOOOOOO 01 

0.100000000 01 
3.100000000 01 

3.100000000 01 
0.988560560 00 
0.973236540 00 
C. 95 756 7060 00 
0.941523710 OO 
0.925083600 OO 
0.908218980 00 
0.990698840 00 
0.873088250 00 
0.054747710 00 
0.935832240 00 
0.016290300 00 

O. 796062360 00 

P. 7750T905D 00 
9. 753258610 00 
0.730504650 00 
0. 706699760 9C 
0.68170I52D 00 
3 . 655332680 00 
0.627368640 00 
0 .59751 t 7 00 00 
0-565389720 00 
0.530442890 00 
0.49.1887190 00 
■“.448*91750 00 
C.39614S96D 00 
0 • 3 366 34 1 3D 00 
rt . 2526747*0 OO 
0.0 

0.252t74T*0 00 
0.316634130 00 
0.398145960 00 
O. 448491 75D 00 
491087190 OO 
O.53O44P09D 00 
0.565309770 00 
9.59751 7 TOD 00 
*5 • 62 736664 D 00 
0.6553326BD 00 
0.681701520 OO 
0,706699760 00 
0.730504650 00 
C. 753250910 00 
0.775079050 00 
0.796062360 00 
0.016290300 00 

n . 835B3224D 00 
0.854747710 00 
0.873088250 00 
0.890898840 00 
0.908218980 00 
0.925063600 00 
0.941523710 00 
0.957567060 00 
0.973238540 00 
0.968560590 00 
0.100000000 01 
0.100000000 01 
0.100000000 01 

o. ioooooooo 01 
O.IOOOOOOOD 01 
O.IOOOOOOOO 01 

o.ioooooooo 01 
o.ioooooooo 01 

0. IOOOOOOOO 01 
O.IOOOOOOOD 01 
0.100000000 01 
0.100000000 01 

O.IOOOOOOOD 01 


Figure B-3. 


Sample output of the WAKE 


p rog ram . 


1 19 


1 


JOUKOtfSK 


CONVERGENCE ATTAINED AFTER 
RESULTS FOR X/C STATION M3 


A ITERATIONS 

X = 0*2*51 5737D 


A IRFOIL/A=0/RN*. 21 /FREE TRANST T I ON/M= . OS/CEOD Y= 


OO 

,03/ 


ETA 

0. 1 0 000 00 OD 01 
0*975000000 00 
0.95000000D 00 
0.92500000D 00 
0.900000000 00 
0*875000000 00 
0*850000000 00 
0*825000000 00 
0*80000000D 00 
0.775000000 00 
0* 750000000 00 
0*725000005 00 
0*700000000 00 
0*875000000 00 
0*650000000 00 
0*625000000 00 
0.600000000 00 
0.575000000 00 
0*550000000 00 
0*525000000 00 
0* 500 0 000 OD 00 
0.47500000D 00 
0* *500 OQOOO 00 
0 • *25000000 00 
0. *00000000 00 
0*375000000 00 
0*350000000 00 
0 .325000000 00 
0*300000000 00 
0*275000000 00 
0*25000 00 OD OO 
0.22500000D 00 
0 * 20000 0000 00 
0.175000000 00 
0*1 50 00 00 OD 00 
0*125000000 OO 
0* L OOOOOOOD 00 
0*750000000-0 I 
O.SOOOOOOOD-Ol 
0*250000000-01 
0*0 

-0* 250000000-01 
-0.500000000-01 
-0. 750000000-01 
-0.100000000 00 
-0*125000000 00 
-0*150000000 00 
-0.17500000D 00 
- 0*200000000 00 
-0*225000000 00 
-0.250000000 00 
-0*275000000 00 

— C* 3 0000 000 O 00 
-0*325000000 00 
-0.350000000 00 
-0*375000000 00 
- 0**00000000 00 
-0*425000000 00 
-0*450000000 00 
-0.475000000 00 
-0*500000000 00 
-0*525000000 00 
-0*550000000 00 
-0.575000000 00 
-0.60000000D 00 
-0*625000000 00 
-0.65000000D 00 
-0*675000000 00 

— 0* 7 OOOOOOOD 00 
-0*725000000 00 
-0*750000000 00 
-0*775000000 00 

— 0* 800 0000 OD 00 
-0*825000000 00 

— 0* 85000 00 OD 00 

— 0*875000000 00 
-0*900000000 00 
-0*925000000 00 

— 0 * 95 000000 D 00 
—0*975000000 00 
- 0.100000000 01 


U/UOEL 

0*1 000000 OD 01 
0*100000000 01 
0.10000000D 01 
0*100000000 01 
0*100000000 01 
0*100000000 01 
0*1 OOOOOOOD 01 
0*100000000 01 
0*1 OOOOOOOD 01 
0*100000000 01 
O.IOOOOOOOD OJ 
0*100000000 01 
0.100000000 01 
0*100000000 01 
0*100000000 01 
0*100000000 01 
0*100000000 01 
0*100000000 01 
0.100000000 01 
0.100000000 01 
O.IOOOOOOOD 0 | 
0*100000000 01 
0.100000000 01 
0*100000000 01 
O.IOOOOOOOD 01 
0*999999960 00 
0.999999630 00 
0.999996970 00 
0.999980160 00 
0*999894480 00 

0*999541050 OQ 
0*998352190 00 
0* 9950 595 SD 00 
0*987462500 00 
0.972668450 00 
0.94807 29 SD 00 
0.91 299 1 800 00 
0-8 70 42 59 2D 00 
0 * 82 788366D 00 
0.7950951 7D 00 
0*783951300 00 
0*7958951 70 00 
0 • 82 788 3860 00 
0*870425920 00 
0*912991800 00 
0.948072950 00 
0 * 97 2666 4 5D 00 
0.987462500 00 
0.995059550 00 
0*990352190 00 
0 * 99954 1 0 5D 00 
0* 9998944 SD 00 
0.999980160 00 
0.999996970 00 
0*999999630 00 
0 *999 9 9996 D 00 
0.100000000 01 
O.IOOOOOOOD 01 
O.IOOOOOOOD 01 
0.100000000 01 
O.IOOOOOOOD 01 
0*100000000 01 
0*100000000 01 
0*100000000 01 
0*100000000 01 
0*100000000 01 
0*100000000 01 
O.IOOOOOOOD 01 
0.100000000 01 
0*100000000 ot 
0*100000000 01 
O.tOOOOOOOO 01 
0*100000000 01 
0*100000000 01 
0.100000000 0| 
0.100000000 01 
0*100000000 01 
0*100000000 ot 
O.IOOOOOOOD 01 
O.IOOOOOOOD 01 
0*100000000 01 


0* 98000 46 00 00 
0*955004600 00 
0* 9 300 04 6 UD 00 
0. 90 500 46 OD 00 
0.860004600 OO 
0*855004600 00 
0.830004600 00 
0*805004600 00 
0.7Q0004600 00 
0 * 7550 046 OD OO 
0*73000*600 00 
0. 705004600 00 
0*680004600 00 
0*655004600 00 
0*630004600 00 
0*605004600 00 
0.580004600 00 
0*555004600 00 
0*530004600 00 
0*505004600 00 
0*480004600 00 
0.455004600 00 
0*430004600 00 
0*405004600 00 
0*380004600 00 
0*355004600 00 
0.330004610 00 
0*305004640 00 
0.28000487D 00 
0.255006170 00 
0.230012280 00 
0 .205035900 00 
0*180111900 00 
0.155318680 00 
0.130799060 00 
0*106767720 OO 
0 . 83*8* I 35D—0 1 
0*61 I62O62B-01 
0*399636620-01 
0*1 9 6 99 5 79 D— 01 
0.0 

—0*1 96995 79D—0 1 
-0. 399636620-01 
-0*61 182062D— 0 1 
-0.634841380-01 
-0.106767720 00 
-0.130799060 00 
-0*155318680 00 

— O* 1 80 111 98D 00 
-0.205035900 00 
-0.230012280 00 
-0*255006170 00 
-0*280004870 00 
-0.305004640 00 
-0.330004610 00 
-0.355004600 00 
-0.38000460D 00 
-0. V 05C 04600 00 
— 0 * 43000460 D 00 

— 0* 4550 0460 D 00 
-0*480004600 00 
-0.505004600 00 
-0*530004600 00 
— 0 *555004600 00 
-0*580004600 00 
— 0* 60 500460D 00 
-0*630004600 00 
-0*655004600 00 
-0*680004600 00 
-0.705004600 00 
-0*730004600 00 
-0*755004600 00 
-0*780004600 00 
-0*805004600 00 
-0*830004600 00 
-0*855004600 00 
-0.880004600 00 
—0 *90 50 0460 D 00 
-0*930004600 00 
-0.955004600 00 
-0*980004600 00 


v/c 

0*1 7326878D 00 
O . I 689 37 060 00 

0*164605340 00 
0.160273620 00 
0.15594I90D 00 
0* 1516101 BD 00 
0*147278 46 D 00 
0*142946740 00 
0. I 386 1 502 0 00 
0 • t 34283300 00 
0*129951580 00 
0*125619860 00 
0*121288140 00 
0.1169S642D 00 
0.1 1 262471 D 00 
0*108292990 00 
0*103961270 00 
0*996295470-01 
0*95 2978270—0 1 
0*909661080-01 
0*866343880-01 
0.823026690-01 
0 * 77970950 D— 01 
0.73639230D— 01 
0.6930751 10-01 
0*6497579 I D— 0 1 
0*606440720-01 
0*563123530-01 
0*51 98063 3D- 0 1 
0*476489 14D-01 
0*433 1 71940-01 
0*389854750-01 
0.346537550-01 
O. 303220 360-01 
0.259903170-01 
0.21 6S8S 970-01 
0. 173268 78D-01 
0.1 2995 1580-01 
0* S6634388D— 02 
0.433171940-02 
0.0 

-0*4331 7194D-02 
-0*866343680-02 
-0. 1 29951 580-01 
-0. 173268700-01 
-0*216585970-01 
-O* 259903170-01 

— O . 3 03220 36D— 0 1 
-0.346537550-01 
-0.389854750-01 
-0.433171940-01 
-0.476489140-01 
-0*519806330-01 
-0*563123530-01 
-0.606440720-01 
-0*649757910-01 
—0-693075 1 10-01 

— 0* 7 36392 30 D— 01 
-0.779709500-0 1 
-0*823026690-01 
-0.666343880-01 
-0.909661080-01 
-0.952978270-01 
-0.996295*70-01 
-0.103961270 00 
-0*108292990 00 
-0.112624710 OO 
-0.1 16956420 00 
-0*121288140 00 
-0.125619860 00 
-0.129951580 00 
— O .134283300 00 
-O. 138615020 00 
-0.142946740 00 
-0*147278460 00 
— 0* 151610180 00 
-0*155941900 00 
—0*1 60 2 7 36 2D 00 
-0* 164605340 00 
-0.168937060 00 
— O • 173266780 00 


«U/UDEL)**2 
O.IOOOOOOOD 01 
O.IOOOOOOOD 01 
0. IOOOOOOOD 01 
0*100000000 Ot 
0.100000000 01 
0.100000000 01 
0*100000000 01 
0-100000000 01 
0.100000000 01 
O.IOOOOOOOD Ot 
O.IOOOOOOOD OJ 
0. I OOOOOOOD 01 
0*100000000 01 
0. t OOOOOOOD 01 
0*100000000 01 
O.IOOOOOOOD 01 
0. 10000000D 01 
0.100000000 01 
0*100000000 01 
0.100000000 01 
o . ioooooood ot 
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ETA 
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0*973000000 00 
0.950000000 00 
C.925COOOOD 00 
0*900000000 00 
0.573000000 00 
0.450000000 00 
0*825000000 00 

o.eoooOooou oo 

0*775000300 00 
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0. 703000000 00 
0*675000000 00 
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0.625000000 00 
0.60000000D 00 
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0.525000000 CO 
0.500000000 00 
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0.350000000 00 
0.32500000D 00 
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0*103000000 00 
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-0.125000000 00 
-0*150000000 00 
-0. 1 75000000 00 
-0*200000000 00 
-0*225000000 00 
-0.250000000 00 
-0.27500000D 00 

- 0. 300000000 00 
-0*325000000 00 
-0.350000COD 00 
-0.375000000 00 
-0*400000000 00 
-0.42500000D 00 
-0*450000000 00 
-0*475000000 00 
-0*500000000 00 
-0.523000000 00 
-0* 550 CO 00 OO 00 
-0*575000000 00 
-0*600000000 00 
—0*62 500 OOOD 00 
-0*650000000 00 
-0*675000000 6 0 

— C* 700000000 00 
-0*72500000D 00 
- 0* 750000000 00 
-0*773000000 00 
-0*800000000 00 
-0*825000000 00 
-0*850000000 00 
—0*873 00 OOOD 00 
-0*900000000 00 
-0*925000000 00 
-0*950000000 00 
-0.975000000 00 
-0*100000000 01 


U/UD6L 

0.100000000 01 
0.100000000 01 
0*100000000 ot 
0*100000000 01 
0*100000000 01 
O.IOCOOOOOO o l 
0.100000000 01 
C.10000000D 01 
0*100000000 01 
0*100000000 01 
O.IOCOOOOOO 01 
O.IOCOOOOOO 01 
0.100000030 01 
0. IOOOOCOOO 01 
0.100000000 01 
O.lQtfOOOOOD Q 1 
0.100000000 01 
0*100000000 01 
0. IOOOOCOOO 0 I 
0*100000000 01 
0.100000000 01 
0*100000000 01 
0. 190000000 01 
0.100000000 01 
C . 1 0 0 0 0 OOOD 01 
0.100000300 01 
0.999999980 00 
0.999999800 00 
0.999998290 00 
0.999988130 00 
0.999932240 00 
0*999681260 00 
0.996756660 00 
0 . 99595 1 5 JO 00 
0* 98892 19BD 00 
0.974360090 00 
0.949511110 OO 
0.914986670 00 
0. 8770 H 1 590 00 
0*846877810 00 
0.835286830 00 
0*846877810 00 
0.877081590 00 
0.914986670 00 
0.949511110 00 
0.974360090 00 
0.988921980 00 
0.99595I5JD 00 
O .996756660 00 
0.999681260 00 
0*999932240 00 
0.9999881 30 00 
0.999998290 00 
0.999999800 00 
0.999999980 00 
0.100000000 01 
0.100000000 01 
0*100000000 01 
0. 1 OOOOOOOD 01 
0*100000000 01 
0*100000000 01 
0.100000000 01 
0*100000000 01 
0* IOOOOCOOO 01 
0*1 OOOOOOOD 01 
0*100000000 01 
0*100000000 01 
0.100000000 01 
0.100000000 01 
0*100000000 01 
O.IOCOOOOOO 01 
0*100000000 01 
0*100000000 01 
0* t OOOOOOOD 01 
0.100000000 01 
0.100000000 01 
0*100000000 01 
0*100000000 01 
0*100000000 01 
0*100000000 01 
0*100000000 01 


0.986592130 00 
0*961592130 00 
0*936592130 00 
0«9L 1592130 00 
C • 386592 l 30 00 
0.861592IJD 00 
C. 836592130 00 
0 .41 1592130 00 
0.786592130 00 
0*761592130 00 
C. 736592130 00 
C.7I 1592 I JO OO 
0.686592130 00 
0.661592130 00 
0.636592130 00 
0.6MS92I3D 00 
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r/c 

0*2474161*40 00 
0*241250110 00 
0.235064230 00 
G.220078J3D 00 
0.222692430 00 
0.216506530 00 
0.210320630 00 
0.204134730 00 
0.19 79408-30 00 
0*191762930 00 
0*1 85577 Cl D 00 
0*179391 I 3D OO 
0.173205220 00 
0*167019320 00 
0*160833420 00 
0. I54647S20 00 

0.148461620 OC 
0.142275720 00 
C. 136069820 00 
0. I 29903920 00 
0.123718020 00 
0.1 1 7532 120 00 
0. I 1 I 340220 00 
0. I 05 1 60 00 
0 .989744 140-01 
0.927885 1 30-01 
C . 8660 ?o ! 2D- 01 
0.804 16* 120—01 
0. 742308 1 10-01 
0.680449 100-01 
0. 618590090-01 
0.556731 080-01 
0.494872070-01 
0.4330 1 3060-01 
0. 371 1 54050-01 
0 .309295040-01 
0*24 74 36,04 0-01 
O. 1 8557 7C 10—0 I 
0. 123718020-01 
0.618590090-02 
0.0 

-0*61 8590 09 D- 02 
-O. 1 2 37 1 BC2D— 0 1 
-0. I 85577030-01 
-0.247436040-01 
-0* 309295040-01 
-0.171 1S4C50-01 
-0.43301 306D-01 
-0 . 494072070-01 
-0.556731080-01 
-0.61 059OC9O-O1 
-0.6804 49 10 0-0 1 
-0.742308 1 I O-Ot 
-0*0041671 20-01 
-0.066026120-01 
-0*927085130-01 
-0*989744140-01 
-0.105160320 00 
-O. 1 1 I 346220 00 
-0*117532120 00 
-0*123718020 00 
-0*129903920 00 
-0*136089820 00 
-0*142275720 00 
-0*148461620 00 
— Q * 15464 752 D 00 

-0*160833420 00 
-0*16701^320 00 
-0*173205220 00 
-0*1 793911 3D 00 
-0*185577030 00 
-0*191762930 00 
-0*197948830 OO 
— 0 *204134 7iD 00 
-0*210320030 00 
-0*216506530 00 
-0*222692430 00 
-0*226878330 00 
-0*233064230 00 
-0*241250130 00 
-0*247436040 00 


Figure B-3. Continued. 


( U/UOEL I 4*2 
o.i oOoooood 01 
0-100000000 01 
o.ioooooooo Ol 
■>. 100000000 01 
0*100000000 ot 
0.10000^500 01 
0.100000000 01 
e.iooccoooo oi 
0.100000000 01 
0.1 OOOOOOOD 6l 
* • lOOCOOOOO 01 

O.IOOOOOOOO 01 

c.iooooooco 01 
0.100000000 01 
0. lOOCOOOOO 01 
C. lOOCOOOOO 01 
O.IOOOOOOOO Ol 
C. lOOCOOOOO 01 
0.100000000 01 
">. 100000000 01 
0. 10000000D 01 

ioooooooo oi 

0.100000000 OI 
0. I OOOOOOOD 01 

lOOCOOOOO Ol 
O.IOOOOOOOO 01 
0.999999960 00 
0.994999600 00 

99999659D 00 
''.999976250 00 
0.999064480 00 
0.999362620 00 
0.997514860 00 
C. 991919440 00 
0,977966600 00 
C.94937758D 00 
0.901571360 00 
0.837200610 00 
0.769272120 00 
0. 717202020 00 
0.697704090 OO 
0.71 7202020 00 
0.769272120 00 
0.837 20061 O 00 
0.901671360 00 
0.94937*390 00 
0.977966680 OC 
0 . 99191 9440 00 
0.997514060 00 
0.999362620 00 
0.999 06 4 4 8 D 00 
0.999976250 00 
0,999996390 00 
0.99999960D 00 
0.999999960 00 
O.IOOOOOOOO 01 
O.IOOOOOOOO 01 
O.IOOOOOOOO 01 
O.IOOOOOOOO 01 
O.IOOOOOOOO 01 
0.100000000 01 
o.ioooooooo ot 
0*100000000 01 
0*100000000 01 
O.IOOOOOOOO 01 
C.IOOCOOOOD Ot 
O.IOOOOOOOO ot 
O.IOOOOOOOO 01 
0*100000000 01 
O.IOOOOOOOO 01 
0*100000000 01 
0.100000000 oi 
0*100000000 01 
O.IOOOOOOOO Ol 
0*100000000 01 
O.IOOOOOOOO 01 
o.ioooooooo 01 
0.100000000 01 
0.1 OOOOOOOD ot 
O.IOOOOOOOO 01 

0*100000000 oi 


CONVERGENCE ATTAINED AFTER 
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3 ITERATIONS 

X = 0-202963*10 
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